restart;
ynp1:=a0*y(tn)+a1*y(tn-h)+a2*y(tn-2*h)+
      h*(bm1*D(y)(tn+h)+b0*D(y)(tn)+b1*D(y)(tn-h)+
         b2*D(y)(tn-2*h))+E5*h^5*(D@@5)(y)(theta)/5!;
r:=y(tn+h)-ynp1;
eq[0]:=eval(subs(y=(x->1),r));
for j from 1 to 5 do
    tmp[j]:=eval(subs(y=(x->x^j),r));
    eq[j]:=coeff(tmp[j],h^j);
    print(eq[j]);
od:
S1:=solve({seq(eq[k],k=0..5)},{a0,bm1,b0,b1,b2,E5});
method:=subs(E5=0,ynp1);
m2:=subs(S1,method);
m3:=eval(subs(D(y)=(x->f(x,y(x))),m2));
f:=(xi,eta)->A*eta;
m4:=y(tn+h)=m3;
ceq:=eval(subs(y=(s->rho^s),m4));
ceq2:=subs({a1=1,a2=1,tn=1,h=1},ceq);
S2:=solve(ceq2,rho);