restart;
# Modified to construct a lower order predictor that is stable...
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))+E4*h^4*(D@@4)(y)(theta)/4!;
r:=y(tn+h)-ynp1;
eq[0]:=eval(subs(y=(x->1),r));
for j from 1 to 4 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]=0,k=0..4)},{a0,b0,b1,b2,E4});
method:=subs(E4=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/3,a1=0,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=-2..1,b=-1..1,contours=[1],
    grid=[100,100],filled=true);
S1;
M1:=subs(S1,a0^2+a1^2+a2^2);
M2:=subs({a1=x,a2=x},M1);
M3:=diff(M2,x);
solve(M3=0,x);
subs({a1=0,a2=1/3},M1); # still less than 1
;
subs({a1=1/3,a2=1/3},M1); # better but also less stable
;