   nh=30;
   Inh=speye(nh);
   T=toeplitz([-2,1,zeros(1,nh-2)]);
   A=kron(T,Inh)+kron(Inh,T);
   I=eye(size(A));
   n=size(A,1);
   C=ones(1,n); C(2:2:end)=-2;
   B=ones(n,1);
   X=care(A,B,C'*C);
   params.m=100;
   params.tol=1e-9;
   params.smin=1e-1; params.smax=1e1;
   params.ch=1;
   params.period=1;
   params.Hritz=1;   % use modified version with spec(Tk - Yk Bk Bk')
   %
   %  X = Z*Z' approx solution to A' X + X A - X B B' X + C' C = 0
   %
   [Z,nrmrestot,VV,K,s]=ARKSM_Riccati(A,B,C,params); 

