ynp1=a0*y[tn]+a1*y[tn-h]+a2*y[tn-2*h]+
h*(bm1*y'[tn+h]+b0*y'[tn]+b1*y'[tn-h]+
b2*y'[tn-2*h])+
E5*h^5*Derivative[5][y][theta]/5!
r=y[tn+h]-ynp1
eq[0]=r==0/.y->Function[x,1]
For[j=1,j<6,j++,
    tmp[j]=r/.y->Function[x,x^j];
    eq[j]=Coefficient[tmp[j],h^j]==0;
    Print[eq[j]]
]
S1=Solve[{eq[0],eq[1],eq[2],eq[3],eq[4],eq[5]},
    {a0,bm1,b0,b1,b2,E5}]
method=ynp1/.E5->0
m2=method/.S1[[1]]
m3=m2/.y'->Function[x,f[x,y[x]]]
f=Function[{xi,eta},A*eta]
m4=y[tn+h]==m3
ceq=m4/.y->Function[s,rho^s]
ceq2=ceq/.{a1->1,a2->1,tn->1,h->1}
S2=Solve[ceq2,rho]
