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