function ps65

[~,P]=lqr([0 1;0 0],[0;1],[1 0;0 0],1,[0;0])
P(1,1)-P(1,2)^2/P(2,2)