function [S,Rts,K]=anasolv(G,cc,dd,x0)
G=ss(G); [Ga,xa]=augment(G,cc,dd,x0);
vec=ones(1,length(Ga.a));
[Aj,T,K,Rts]=jordan_m(Ga.a,vec); 
B=inv(T)*xa;  F=[]; G0=[]; j=1; 
for i=1:length(K), 
   if (imag(Rts(j))>1e-5), 
      F1=[B(j),B(j+1); B(j+1),-B(j)]; 
      G1=diag([1,1]); j=j+2;
   else
      Xx=B(j:j+K(i)-1)'; 
      F1=Xx; G2=1;
      for k=2:K(i)
         G2(k)=G2(k-1)/(k-1);
         Xx=[Xx(2:K(i)), 0];
         F1=[F1; Xx];
      end
      G1=diag(G2); j=j+K(i); 
   end
   mm=size(F1); nn=length(F); 
   F=[F,zeros(nn,mm(2)); zeros(mm(1),nn),F1];
   G0=[G0, zeros(nn,mm(2)); 
      zeros(mm(1),nn), G1];
end
S=T*F*G0;

%jordan_m function is the supporting function for 
%anasolv().  If performs Jordan realization to
%the original system model.
function [AJord,T,JJ,Rts]=jordan_m(A,V0,Eps)
if nargin==2, Eps=1e-5; end
Rts=eig(A); Rts=sort(Rts);
Rts=Rts(length(Rts):-1:1); V=[]; 
Ndim=length(Rts); M=V0;
for i=2:Ndim, M=[M; M(i-1,:)*A]; end
k = 0; JJ = [];
while k < Ndim,
   k=k+1; K=0; key=1; 
   Rts0=Rts(k); Rts1=Rts(k);
   if (abs(imag(Rts0))>Eps), 
      key=2; K=K+1; k=k+2;
   else
      while (k<=Ndim & abs(Rts0-Rts1)<1e-5),
        k=k+1; K=K+1; key=0;
        if k<=Ndim, Rts1=Rts(k); end   
   end, end
   if key~=1, JJ=[JJ K]; k=k-1; end
end
V=[]; K=0;
for i=1:length(JJ),  
   V0=[]; Rts0=Rts(K+1);
   if (abs(imag(Rts0))<eps),  
      K=K+JJ(i);
      for j=1:Ndim,
         V0=[V0; real(Rts0)^(j-1)];
      end
      if JJ(i)>1,
         V1=V0; 
         for k=2:JJ(i), 
            V2=zeros(Ndim,1);
            for j=k:Ndim
               V2(j)=V1(j-1)*(j-1)/(k-1);
            end
            V0=[V0 V2]; V1=V2;
      end, end
   else, K=K+2;
      for j=1:Ndim,
         V0=[V0; 
           real(Rts0^(j-1)) imag(Rts0^(j-1))];
   end, end
   V=[V V0];
end
T=inv(M)*V; AJord=inv(T)*A*T;
for i=1:Ndim, for j=1:Ndim
   if abs(AJord(i,j))<Eps, AJord(i,j)=0; end
end, end