restart; cor1 := (1 - a1 - a2)*y(tn) + a1*y(tn - h) + a2*y(tn - 2*h) + h*((3/8 - a1/24)*D(y)(tn + h) + (19/24 + (13*a1)/24 + a2/3)*D(y)(tn) + (-5/24 + (13*a1)/24 + (4*a2)/3)*D(y)(tn - h) + (1/24 - a1/24 + a2/3)*D(y)(tn - 2*h)); pred1 := (-8 - a2)*y(tn) + 9*y(tn - h) + a2*y(tn - 2*h) + h*((17/3 + a2/3)*D(y)(tn) + ((4*a2)/3 + 14/3)*D(y)(tn - h) + (-1/3 + a2/3)*D(y)(tn - 2*h)); c2:=subs({a1=0,a2=2/5},cor1); p2:=subs(a2=1,pred1); c3:=eval(subs(D(y)=(x->f(x,y(x))),c2)); p3:=eval(subs(D(y)=(x->f(x,y(x))),p2)); pc3:=subs(y(tn+h)=p3,c3); f:=(xi,eta)->A*eta; m4:=y(tn+h)=pc3; ceq:=eval(subs(y=(s->rho^s),m4)); ceq2:=subs({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=-0.75..0.25,b=-0.5..0.5,contours=[1], grid=[100,100],filled=true);