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];
I=eye(4);O=0*I; B=I; Q=I; R=I; N=10000; L=zeros(4)
%
J=[];
for i=1:20;
x=zeros(4,1); w=[]; p=zeros(4,1); 
for t=1:N, w=[w randn(4,1)]; x=[x (A-B*L)*x(:,t)+w(:,t)]; end,
for t=1:N, p=[(A-B*L)'*p(:,1)-2*Q*x(:,N+2-t) p]; end,
xx=x*x'/N;
ww=w*w'/N;
px=p*x'/N;
J=[J trace(xx+L*xx*L')];
Ji=J(i)
DL=2*R*L*xx+B'*px;
DL(1,3)=0;DL(1,4)=0;DL(2,4)=0;DL(3,1)=0;DL(4,1)=0;DL(4,2)=0;
L=L-.02*DL
end,
plot(J(1:20),'o')