G1=tf([5,100],[1,8,32,80,100]); G=ss(G1);
u0=[2,-1,-24]; Y0=[0.3; 0.5; 0; 0];
hatY0=Y0-[0; G.c*G.b*u0(1); 
     G.c*G.b*u0(2)+G.c*G.a*G.b*u0(1); 
     G.c*G.b*u0(3)+G.c*G.a*G.b*u0(2)+G.c*G.a^2*G.b*u0(1)];
Wo=obsv(G); X0=inv(Wo)*hatY0;  
[Ga,Xa]=augment(G,1,[],X0);

T=0.05; AA=expm(Ga.a*T); t=[0]; y=[Ga.c*Xa]; x=Xa;
for i=1:100
   x=AA*x; y=[y; Ga.c*x]; t=[t; i*T];
end
plot(t,y); axis([0,5,0,1.2])