clf;
A=rand(2,2);B=rand(2,1);   //two states, one input
Q=diag([2,5]);R=2;     //Usual notations x'Qx + u'Ru
Big=sysdiag(Q,R);    //Now we calculate C1 and D12
[w,wp]=fullrf(Big);C1=wp(:,1:2);D12=wp(:,3:$);   //[C1,D12]'*[C1,D12]=Big
P=syslin('c',A,B,C1,D12);    //The plant (continuous-time)
[K,X]=lqr(P);
K
spec(A+B*K) ;   //check stability
K = 10000*K;
U = A+B*K;
P2 = syslin('c',U,B,C1,D12);
X0 =[1;1];
t = 0:0.01:100;
u = zeros(t);
[y,x] = csim(u,t,P2,X0);
plot(t,x(1,:),t,x(2,:));
halt;
clf;
plot(x(1,:),x(2,:));