restart;
ynp1:=a0*y(tn)+a1*y(tn-h)+a2*y(tn-2*h)+
      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,a1,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({a2=1,tn=1,h=1},ceq);
S2:=solve(ceq2,rho):
Z1:=subs(A=a+I*b,abs(S2[1])):
Z2:=subs(A=a+I*b,abs(S2[2])):
Z3:=subs(A=a+I*b,abs(S2[3])):
with(plots):
contourplot(max(Z1,Z2,Z3),a=-3..3,b=-2..2,contours=[5,6,7,8,9,10],
    grid=[100,100],filled=true);
# Note that max is >4 for all values of a and b.
# Values of blue >10
pred1:=subs(S1,subs(E5=0,ynp1));