%PLANT
al = [1+.3*s, .5*s; .1*s 1+0.7*s];
bl = [.2 .4;.6 .8];
bl = pol(bl);
[a,b,c,d] = lmf2ss(bl,al);
G = ss(a,b,c,d);

%CONTROLLER
als = al*s;
[As,Bs,Cs,Ds] = lmf2ss(bl,als);
    [Brs,Ar] = ss2rmf(As,Bs,eye(4));
Q =10*eye(4);
Fr = spf(Ar'*Ar+Brs'*Q*Brs);
[Xl,Yl] = xaybc(Ar,Brs,Fr);
K=mldivide(Xl,Yl);
K=lcoef(K);
[Bls,Al] = ss2lmf(As,eye(4),Cs);
%V=10*eye(4);
%Ol = spf(Al*Al'+Bls*V*Bls');
Ol =[-3.1-3.4*s-0.18*s^2    -0.5-0.58*s-0.011*s^2;
     -0.5-0.62*s-0.055*s^2   3.1+4*s+0.9*s^2]; 
 
[Xr,Yr] = axbyc(Al,Bls,Ol);
L=mrdivide(Yr,Xr);
L=lcoef(L);
acp = As-L*Cs-Bs*K;
bcp = L;
ccp = K;
dcp = zeros(2);

[brp,arp]=ss2rmf(acp,bcp,ccp,dcp);
[ac,bc,cc,dc] = rmf2ss(brp,arp*s)