function [X,Xd,Xdd,t] = cendiff(M,C,K,f,dt,Xi,Xdi)  
%================================================================  
%function [X,t] = newmark(M,C,K,f,dt,gamma,beta,Xi,Xdi)  
%  
%  (c) Ajax, (all rights reserved)  
%  Central South University  
%  October 13, 2006  
%  
%  This function is intended to perform the numerical  
%  integration of a structural system subjected to an  
%  external dynamic excitation such as a wind or earthquake.  
%  The structural model is assumed to be a lumped mass shear  
%  model.The integration scheme utilized for this analysis  
%  is the central difference method. The central difference  
%  method is an implicit time steping scheme so stability of  
%  the system need not be considered. 
%  
%  Input Variables:  
%  [M] = Mass Matrix (nxn)  
%  [C] = Damping Matrix (nxn)  
%  [K] = Stiffness Matrix (nxn)  
%  {f} = Excitation Vector (mx1)  
%  dt = Time Stepping Increment  
%  Xi = Initial Displacement Vector (nx1)  
%  Xdi = Initial Velocity Vector (nx1)  
%  
%  Output Variables:  
%  {t} = Time Vector (mx1)  
%  [X] = Response Matrix (mxn)  
%================================================================  
 
  n = size(M,1);  
  fdimc = size(f,2);  
  fdimr = size(f,1);  
  X(:,2)=Xi; 
  Xd(:,1)=Xdi; 
   
%  Check Input Excitation  
%  ======================  
   if(fdimc==n) 
       f=f'; 
   end 
   m=size(f,2); 
    
%  Coefficients  
%  ============  
        c0 = 1/(dt*dt) ;  
        c1 = 1/(2*dt) ;  
        c2 = 2*c0 ;  
         
%Initialize Mass Matricies  
% ==============================  
        Meff = c0*M + c1*C ;  
        Minv = inv(Meff) ; 
         
% Initial Acceleration 
% ====================  
        R0=f(:,1); 
        Xdd(:,1)= inv(M)*(f - K*Xi - C*Xdi); 
 
% Perform First Step  
% ================== 
        X(:,1) = X(:,2) - 1/(2*c1)*Xd(:,1) + 1/(2*c0)*Xdd(:,1) ;  
        f(:,1) = f(:,1) - (K - c2*M)*X(:,2) - (c0*M - c1*C)*X(:,1); 
        X(:,3) = Minv*f(:,1); 
        Xd(:,2) = c1*(X(:,3) - X(:,1)); 
        Xdd(:,2) = c0*(X(:,3) - 2*X(:,2) + X(:,1)); 
         
%Perform Subsequent Steps  
% ========================  
        for i=1:size(f,2)-1;  
         X(:,i+1) = X(:,i+2) - 1/(2*c1)*Xd(:,i+1) + 1/(2*c0)*Xdd(i+1) ;  
         f(:,i+1) = f(:,i+1) - (K - c2*M)*X(:,i+2) - (c0*M - c1*C)*X(:,i+1); 
         X(:,i+3) = Minv*f(:,i+1); 
         Xd(:,i+2) = c1*(X(:,i+3) - X(:,i+1)); 
         Xdd(:,i+2) = c0*(X(:,i+3) - 2*X(:,i+2) + X(:,i+1)); 
        end;  
 
          
%  Strip Off Padded Response Zeros  
%  ===============================  
  X = X' 
  Xd=Xd' 
  Xdd=Xdd' 
   
%  Generate the Time Vector  
%  ========================  
  for i=0:1:size(X,1)-1  
    t(i+1,1) = i*dt;  
  end;  
  t