A=[0 1 0 0; 0 0 1 0; 0 0 0 1; 10e-11 -24010 -2539 -54];
B=[0; 0; 0; 31940];
C=[1 0 0 0];
D=[0];

Td=0.01;

sys=ss(A,B,C,D);
dsys=c2d(sys,Td,'zoh');
[M,N,C,D]=ssdata(dsys);

Q=0.01*eye(4);
R=0.05;

N=50;
x=zeros(4,N);
Po=randn*eye(4,4);
y=zeros(1,N);
u=zeros(1,N);
e=randn(N);
e=e(1,:);

S=zeros(4)

x(1)=0;
P(1)=0;
i=2;
for i=2:N
   L=Po*C'*inv(R+C*Po*C');
%  y(i)=C*x(i)+D*u(i)+e(i);
   x(:,i)=(x(i-1)+L(i)*(y(i)-C*x(i-1)-D*u(i)))';
   Po=Po*A'+Q-(A*Po*C')*inv(C*Po*C'+R)*(C*Po*A');
end;