(* ab2stab.wl -- The linear stability domain for an AB2 *)

method=y[t+2*h]==y[t+h]+h*(3/2*f[t+h,y[t+h]]
-1/2*f[t,y[t]])
linear=method/.f->Function[{t,y},lambda*y]
eq1=Simplify[linear/.{h->1,lambda->z}]
eq2=eq1/.y->Function[t,w^t]
eq3=eq2/.t->0
roots=Solve[eq3,w]
z1=Abs[w]/.roots[[1]]
z2=Abs[w]/.roots[[2]]
z1ab=z1/.{z->a+I*b}
z2ab=z2/.{z->a+I*b}
ContourPlot[Max[z1ab,z2ab],{a,-2.5,0.5},{b,-1.5,1.5},
Contours->{1},ContourShading->{Red,Blue}]
