RK:=proc(s,p)
  local TaylorPhi,RungeKuttaPhi,d,vars,eqns,k,i,j,dt;
  global a,b,c;
   D(y):=t->f(t,y(t)):                      # Taylor series
  TaylorPhi:=convert(taylor(y(t+dt),dt=0,p+1),polynom):
  TaylorPhi:=normal((TaylorPhi-y(t))/dt);
  c[1]:=0:                                  # RK-ansatz
  for i from 1 to s do
   k[i]:=taylor(f(t+c[i]*dt,y(t)+sum(a[i,j]*k[j],j=1..i-1)*dt),dt=0,p):
  od:
  RungeKuttaPhi:=sum(b[j]*k[j],j=1..s):
  RungeKuttaPhi:=series (RungeKuttaPhi,dt,p):
  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)};
                                            # simplifying assumptions:
  eqns:=eqns union {seq(sum(a[i,'j'],'j'=1..i-1)-c[i],i=2..s)}; 
  solve(eqns,vars);                                          
end;