restart;
k[1]:=f(y(t0)):
k[2]:=taylor(f(y(t0)+dt*a[2,1]*k[1]),dt,5):
k[3]:=taylor(f(y(t0)+dt*(a[3,1]*k[1]+a[3,2]*k[2])),dt,5):
k[4]:=taylor(f(y(t0)+dt*(a[4,1]*k[1]+a[4,2]*k[2]+a[4,3]*k[3])),dt,5):

y1:=taylor(y(t0)+dt*sum(b[j]*k[j],j=1..4),dt,5):
tau:=taylor(y(t0+dt),dt,5)-y1:
D(y):=t->f(y(t));
alias(f=f(y(t0)),f[y]=D(f)(y(t0)),f[yy]=(D@@2)(f)(y(t0)),f[yyy]=(D@@3)(f)(y(t0)));
p:=collect(convert(tau,polynom),dt);
o1:=collect(op(4,p)/dt/f,[f,f[y],f[yy],f[yyy]]);
o2:=collect(op(3,p)/dt^2/(f[y]*f),[f,f[y],f[yy],f[yyy]]);
o3:=collect(op(2,p)/dt^3,[f,f[y],f[yy],f[yyy]],factor);
collect(op(1,o3)/(f[yy]*f^2),[b[1],b[2],b[3],b[4]],factor);
collect(op(2,o3)/(f[y]^2*f),[b[1],b[2],b[3],b[4]],factor);

o4:=collect(op(1,p)/dt^4,[f,f[y],f[yy],f[yyy]]);
collect(op(1,o4)/(f[yyy]*f^3),[b[1],b[2],b[3],b[4]],factor);
collect(op(2,o4)/(f[yy]*f^2*f[y]),[b[1],b[2],b[3],b[4]],factor);