Euler approximation expansion

Clear[p, eq1, eq2]

p[x_] := (Subscript[a, 0]) + (Subscript[a, 1]) (x - b)

feuler[x_, b_] := f[b] + (x - b) f'[b]

eq1 = p[b] == feuler[b, b]

eq2 = p'[b] == D[feuler[x, b], x] /. x -> b

Solve[{eq1, eq2}, {(Subscript[a, 0]), (Subscript[a, 1])}]

Trapezoid approximation expansion

Clear[p, eq1, eq2, eq3]

p[x_] := 
      (Subscript[a, 0]) + (Subscript[a, 1]) (x - 
              b) + (Subscript[a, 2]) (((x - b))^2) + (Subscript[a, 3]) (((x - b))^3)

ftrap[x_, b_] := 
      f[b] + (1/2) (x - b) (f'[b] + f'[x])

eq1 = p[b] == ftrap[b, b]

eq2 = p'[b] == D[ftrap[x, b], x] /. x -> b

eq3 = p''[b] == D[ftrap[x, b], {x, 2}] /. 
        x -> b

eq4 = p'''[b] == D[ftrap[x, b], {x, 3}] /. 
        x -> b

Solve[{eq1, eq2, eq3, eq4}, {(Subscript[a, 0]), (Subscript[a, 1]), (Subscript[a, 2]), (Subscript[a, 3])}]

Midpoint approximation expansion

Clear[p, eq1, eq2, eq3, eq4]

p[x_] := 
      (Subscript[a, 0]) + (Subscript[a, 1]) (x - 
              b) + (Subscript[a, 2]) (((x - b))^2) + (Subscript[a, 3]) (((x - b))^3)

fmidpt[x_, b_] := f[b] + (x - b) f'[(((b + x))/2)]

eq1 = p[b] == fmidpt[b, b]

eq2 = p'[b] == D[fmidpt[x, b], x] /. x -> b

eq3 = p''[b] == D[fmidpt[x, b], {x, 2}] /. 
        x -> b

eq4 = p'''[b] == D[fmidpt[x, b], {x, 3}] /. 
        x -> b

Solve[{eq1, eq2, eq3, eq4}, {(Subscript[a, 0]), (Subscript[a, 1]), (Subscript[a, 2]), (Subscript[a, 3])}]

Runge-Kutta approximation expansion

Clear[p, eq1, eq2, eq3, eq4, eq5, eq6]

p[x_] := 
      (Subscript[a, 0]) + (Subscript[a, 1]) (x - 
              b) + (Subscript[a, 2]) (((x - b))^2) + (Subscript[a, 3]) (((x - b))^3) + (Subscript[a, 4]) (((x - b))^4) + (Subscript[a, 5]) (((x - b))^5)

frunge[x_, 
        b_] := (1/3) (2  fmidpt[x, b] + ftrap[x, b])

eq1 = p[b] == frunge[b, b]

eq2 = p'[b] == D[frunge[x, b], x] /. x -> b

eq3 = p''[b] == D[frunge[x, b], {x, 2}] /. 
        x -> b

eq4 = p'''[b] == D[frunge[x, b], {x, 3}] /. 
        x -> b

eq5 = p''''[b] == D[frunge[x, b], {x, 4}] /. 
        x -> b

eq6 = p'''''[b] == D[frunge[x, b], {x, 5}] /. 
        x -> b

Solve[{eq1, eq2, eq3, eq4, eq5, eq6}, {(Subscript[a, 0]), (Subscript[a, 1]), (Subscript[a, 2]), (Subscript[a, 3]), (Subscript[a, 4]), 
        (Subscript[a, 5])}]