Oo=zeros(3);
Q=[-0.4 0.2 0.2; 0.2 (-0.5) 0.3; 0.6 0.3 (-0.9)];
D=[0 0.1 0.2; 0.2 0 0.1;0.1 0.2 0];
L=[0.1 0 0; 0 0.2 0; 0 0 0.3];
mu=1;
sigma=0.01;
N=40;
I=eye(3);
S=[0.1 0.02 0.04; 0.04 0.2 0.03; 0.06 0.06 0.3];
A=zeros(6);
A(4:6,4:6)=S;
B=zeros(6);
B(1:3,4:6)=S;
B(4:6,1:3)=mu*I;
B(4:6,4:6)=Q-S-mu*I;
V=zeros(6);

G=zeros(246);
V(1:3,4:6)=0*sigma*I;
B(1:3,1:3)=Q-S-0*sigma*I;
G(1:6,1:6)=B;
G(1:6,7:12)=A;
j=1;
for i=7:6:(N-1)*6+1
    V(1:3,4:6)=j*sigma*I;
    B(1:3,1:3)=Q-S-j*sigma*I;
    G(i:(i+5),(i-6):(i-1))=V;
    G(i:(i+5),i:(i+5))=B;
    G(i:(i+5),(i+6):(i+11))=A;
    j=j+1;
end;
V(1:3,4:6)=N*sigma*I;
B(1:3,1:3)=Q-S-N*sigma*I;
G(241:246,235:240)=V;
G(241:246,241:246)=B;

M=zeros(246,247);
M(1:246,1:246)=G;
M(1:246,247)=1;

bb=zeros(1,247);
bb(1,247)=1;

p=bb*M'/(M*M');
p1=zeros(1,40);
for i=1:1:40
    p1(i)=p((i*6)-5);
end;
plot(p1);