A=[   0.6 0.1 0   0
      0.3 0.6 0.1 0
      0   0.3 0.6 0.1
      0   0   0.3 0.6];
Ad=diag([.6 .6 .6 .6]);
At=A-Ad;
I=eye(4);O=0*I;
B=I;
Q=I;
R=I;
[P,poles,L]=dare(A,B,Q,R);
M=P*(A-B*L);
X=dlyap(A-B*L,I);
J=trace(X+L*X*L');
Z=X+L*X*L'-2*M*X*At'+2*At'*M*X;
dZ=diag(Z)
sdZ=sum(dZ)
%Lb=L; Lb(1,3)=0;Lb(1,4)=0;Lb(2,4)=0;Lb(3,1)=0;Lb(4,1)=0;Lb(4,2)=0
Lb=diag(diag(L));
Mb=Lb;
Xb=dlyap(A-B*Lb,I);
Jb=trace(Xb+Lb*Xb*Lb');
Zb=Xb+Lb*Xb*Lb'-2*Mb*Xb*At'+2*At'*Mb*Xb;
dZb=diag(Zb)
sdZb=sum(dZb)
%
Ae=mdiag(A-B*Lb,Ad);
Be=[O O; B I];
Qe=[O Mb'*At';At*Mb Q];
Re=mdiag(R,O);
Se=[O -Mb';O O];
[Pe,ei,Le]=dare(Ae,Be,Qe,Re,Se);
Xe=dlyap(Ae-Be*Le,[I I;I I]);
Zl=[O I]*Xe*[O I]'+Le(1:4,:)*Xe*Le(1:4,:)'+2*[Mb O]*Xe*Le(5:8,:)'+2*At'*[Mb O]*Xe*[O;I];
dZl=diag(Zl)
disp('Warning: This script works in Matlab 6, but seems to give wrong numbers in Matlab 8')
