# Find the Linear Domain of Stability for an RK Method
restart;
with(LinearAlgebra):
A:=Matrix([[0,0],[1/3,1/3]]);

b:=Vector([1/4,3/4]);

n:=Dimension(b);
Id:=Matrix(n,n,shape=identity);
T1:=Id-z*A;
T2:=MatrixInverse(T1);
onesvector:=Vector(n,1);
T3:=Multiply(T2,onesvector);
T4:=1+z*Multiply(Transpose(b),T3);
r:=simplify(T4);








with(plots):
rxy:=subs(z=x+I*y,r);
contourplot(abs(rxy),x=-10..10,y=-10..10,contours=[1],
   numpoints=1000,axes=boxed,view=[-10..10,-10..10]);