RungeKutta := proc(s, m) local TaylorPhi, RungeKuttaPhi, d, vars, eqns, k, i, j; global a, b, c, h; # Taylor series D(y) := x -> f(x,y(x)): TaylorPhi := convert(taylor(y(x+h),h=0,m+1), polynom): TaylorPhi := normal((TaylorPhi - y(x))/h); # RK-Ansatz: c[1] := 0; for i from 1 to s do k[i] := taylor(f(x+c[i]*h, y(x) + sum(a[i,j]*k[j], j=1..i-1)*h), h=0, m): od: RungeKuttaPhi := 0: for i from 1 to s do RungeKuttaPhi := RungeKuttaPhi + b[i] * k[i]: od: RungeKuttaPhi := series (RungeKuttaPhi, h, m): RungeKuttaPhi := convert(RungeKuttaPhi, polynom); d := expand(TaylorPhi - RungeKuttaPhi): vars := {seq(c[i], i=2..s), seq(b[i], i=1..s), seq((seq(a[i,j], j=1..i-1)), i = 2..s)}; eqns := {coeffs(d, indets(d) minus vars)}; # symmetry condition: eqns := eqns union {seq(sum(a[i,'j'], 'j'=1..i-1) - c[i], i=2..s)}; solve(eqns, vars); end: