restart;
yp:=y(xn-h)+2*h*f(xn,y(xn));
ynp1:=y(xn)+h/2*(f(xn,y(xn))+f(xn+h,yp));
f:=(xi,eta)->A*eta;
method:=y(xn+h)=ynp1;
ceq:=eval(subs(y=(s->rho^s),method));
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[1]));
Z2:=subs(A=a+I*b,abs(S[2]));
with(plots):
contourplot(max(Z1,Z2),a=-3..1,b=-2..2,contours=[1],grid=[100,100],filled=true);