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}]