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]