(* irkgauss.wl -- The linear stability of Gauss-Legendre IRK *) d1=Sqrt[3]/6 eq1=k1==f[t+(1/2-d1)*h,y[t]+1/4*h*k1+(1/4-d1)*k2] eq2=k2==f[t+(1/2+d1)*h,y[t]+(1/4+d1)*h*k1+1/4*h*k2] method=y[t+h]==y[t]+h*(1/2*k1+1/2*k2) f=Function[{t,y},lambda*y] eq1m=eq1/.{h->1,lambda->z} eq2m=eq2/.{h->1,lambda->z} methodm=method/.{h->1,lambda->z} y=Function[t,w^t] eq1w=eq1m/.t->0 eq2w=eq2m/.t->0 methodw=methodm/.t->0 s1=Solve[{eq1w,eq2w,methodw},{k1,k2,w}] wabs=(Abs[w]/.s1)[[1]] w1abs=wabs/.z->a+I*b ContourPlot[w1abs,{a,-1.5,1.5},{b,-1.5,1.5}, Contours->{1},ContourShading->{Red,Blue}]