(* Question 4 *)
ynp1=y[tn]+h/2*(f[tn,y[tn]]+f[tn+h,y[tn+h]])
f=Function[{t,y},A*y]
method=y[tn+h]==ynp1

ceq=method/.y->Function[s,rho^s]
ceq2=ceq/.{tn->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}

(* the red area is linearly stable, the blue unstable *)ContourPlot[Z1,{a,-2,2},{b,-2,2},Contours->{1},ContourShading->{Red,Blue},PlotRange->All,Exclusions->None]