A=[0 1 3 -1;0 -1 -2 2;-1 0 -2 -1;-1 0 -3 0]
b1=[1;0;0;0]
b2=[0;0;1;0]
c1=[0 1 0 0]
c2=[1 0 0 0]
eig(A)
b=b1
c=c2
rank([b,A*b,A*A*b,A*A*A*b])
rank([c;c*A;c*A*A;c*A*A*A])
H=ss(A,b,c,0);
tf(H)
syms s
factor(c*inv(s*eye(4)-A)*b)