A = [0 1 0 0; 0 0 1 0; 0 0 0 1; -10^(-11) -24010 -2539 -54];
B=[0; 0; 0; 31940];
C=eye(4);
D=[0];
T=0.01;

disk = ss(A,B,C,D,T);

Q=[1000 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1];
R=1000;

P=Q;
N=1000;
for k=N:-1:0
   P=A'*P*A + Q - A'*P*B*inv(R+B'*P*B)*B'*P*A;
end

K=inv(R+B'*P*B)*B'*P*A
eig(A-B*K)