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