% Exercise 12 of the Elgersburg School 2010
% Volker Mehrmann and Stephan Trenn
clear all;
syms R L C RL RLp RLpp RLppp f fp fpp
%x=(iR,vR,iC,vC,iL,vL,iRL,vRL)
A=sym([R,-1,0,0,0,0,0,0;0,0,1,0,0,0,0,0;0,0,0,0,0,1,0,0;0,0,0,0,0,0,RL,-1;1,0,0,0,0,0,0,0;1,0,-1,0,-1,0,0,0;0,0,0,0,1,0,-1,0;0,1,0,1,0,0,0,0;0,0,0,-1,0,1,0,1;0,0,0,0,0,0,0,1]);
B=sym([0,0;0,0;0,0;0,0;-1,0;0,0;0,0;0,1;0,0;0,0]);
cA=[A,B];
cAp=sym(zeros(10,10));cAp(4,7)=RLp;
cApp=sym(zeros(10,10));cApp(4,7)=RLpp;
cAppp=sym(zeros(10,10));cAppp(4,7)=RLppp;
E=sym(zeros(10,8)); E(2,4)=C;E(3,5)=L;
cE=[E,zeros(10,2)];
M0=cE;
N0=cA;
Z230=null(M0');
rank(Z230'*N0)
size(Z230'*N0)
M1=[cE,zeros(10,10);-cA,cE];
N1=[cA,zeros(10,10);cAp,zeros(10,10)];
Z231=null(M1');
rank(Z231'*N1)
size(Z231'*N1)
M2=[cE,zeros(10,20);-cA,cE,zeros(10,10);-cAp,-cA,cE];
N2=[cA,zeros(10,20);cAp,zeros(10,20);cApp,zeros(10,20)];
Z232=null(M2');
rank(Z232'*N2)
size(Z232'*N2)
M3=[cE,zeros(10,30);-cA,cE,zeros(10,20);-cAp,-cA,cE,zeros(10,10);-cApp,-cAp,-cA,cE];
N3=[cA,zeros(10,30);cAp,zeros(10,30);cApp,zeros(10,30);cAppp,zeros(10,30)];
Z233=null(M3');
rank(Z233'*N3)
size(Z233'*N3)

Z2=Z232;
A2hat=Z2'*N2*[eye(10,10);zeros(20,10)];
A22=A2hat(1:8,1:8);
B2=A2hat(1:8,9:10);
Zhat=null(transpose(A2hat(:,1:8)));
B3=transpose(Zhat)*A2hat*[zeros(8,2);eye(2,2)];

%doesn't work for some reason:
%B3^(-1)*transpose(Zhat)*transpose(Z23)*[0;0;0;0;0;0;0;0;0;f;0;0;0;0;0;0;0;0;0;fp;0;0;0;0;0;0;0;0;0;fpp];

H=simplify(B3^(-1)*transpose(Zhat)*transpose(Z2));

%also doesn't work:
%f*H(:,10)+fp*H(:,20)+f*H(:,30)

H(:,10) %coefficient of f
H(:,20) %coefficient of f'
H(:,30) %coefficient of f''