yp=y[xn-h]+2*h*f[xn,y[xn]]
ynp1=y[xn]+h/2*(f[xn,y[xn]]+f[xn+h,yp])
f=Function[{xi,eta},A*eta]
method=y[xn+h]==ynp1
ceq=method/.y->Function[s,rho^s]
ceq2=ceq/.{xn->0,h->1}
S=Solve[ceq2,rho]
# the linear stability region is all values of A such that |rho|<1
Z1t=Abs[rho]/.S[[1]]
Z1 = Z1t/.{A->a+I*b}
Z2t=Abs[rho]/.S[[2]]
Z2=Z2t/.{A->a+I*b}
ContourPlot[Max[Z1,Z2],{a,-3,1},{b,-2,2},Contours->{1},ContourShading->{Red,Blue}]
