restart;
eq:=D(y)=(s->f(s,y(s)));
yp:=y(xn)+h*f(xn,y(xn));
ynp1:=y(xn)+h/2*(f(xn,y(xn))+f(xn+h,yp));
r:=y(xn+h)-ynp1;
subs(h=0,r);
T1:=diff(r,h);
dr:=eval(subs(eq,T1));
subs(h=0,dr);
ddr:=eval(subs(eq,diff(dr,h)));
subs(h=0,ddr);
d3r:=eval(subs(eq,diff(ddr,h))):
simplify(subs(h=0,d3r));
d4r:=simplify(eval(subs(eq,diff(d3r,h)))):
simplify(subs(h=0,d4r));
f:=(xi,eta)->A*eta;
method:=y(xn+h)=ynp1;
ceq:=eval(subs(y=(s->rho^s),method));
solve(ceq,rho);
ceq2:=subs({xn=0,h=1},ceq);
S:=solve(ceq2,rho);
# the linear stability region is all values of A such that |rho|<1
Z1:=subs(A=a+I*b,abs(S));
with(plots):
contourplot(Z1,a=-2..2,b=-2..2,contours=[1,.9,.8],grid=[100,100]);
contourplot(Z1,a=-3..1,b=-2..2,contours=[1],grid=[100,100],filled=true);
R:=S/exp(A);