nh=30;
T=diag(2*ones(nh,1))+diag(-ones(nh-1,1),1)+diag(-ones(nh-1,1),-1)+diag(ones(nh-2,1),-2);
I=speye(nh);
A=-(kron(T,I)+kron(I,T));
n=nh^2;
B=randn(n,2);
m=100;
tol=1e-9;
tolY=1e-12;
opts.tol=1e-2;
s1=eigs(-A,1,'largestreal',opts)
s2=eigs(-A,1,'smallestreal',opts)
params.tol=1e-8;
params.smax=s1;
params.smin=s2;
params.m=100;
params.ch=1;
params.period=1;
Z=RKSM_Lyap_real(A,B,params);

fprintf('final true absolute residual norm: \n')
disp(norm(A*Z*Z'+Z*Z'*A'+B*B'))    %this matrix should never be formed for n large 
