% CENTRAL DIFFERENCE
% femem_springdyn_cd3
% damping - 2DOF
clear
   nstep=100;
   T=10;
   %T0=0.2;
   F0=1000;
   M=[1 0;0 .1];
   k=10;
   k0=100;
   K=[k0+k -k;-k k];
    
   dt=T/nstep;
   t=[1:nstep]*dt;

   
   NumEQ=2;
 
   F=zeros(NumEQ,nstep);
   V=zeros(NumEQ,nstep);
   U=zeros(NumEQ,nstep);
   AC=zeros(NumEQ,nstep);

   
   w=10;
   F(1,:)=F0*sin(w*t);
  
   
   AC(:,1)=M\(F(:,1)-K*U(:,1));
   UP=U(:,1)-dt*V(:,1)+dt^2/2*AC(:,1);
   U(:,2)=M\(dt^2*F(:,1)+(2*M-dt^2*K)*U(:,1)-M*UP);
   
  for n=3:nstep
      U(:,n)=M\(dt^2*F(:,n-1)+(2*M-dt^2*K)*U(:,n-1)-M*U(:,n-2));
      AC(:,n-1)=M\(F(:,n-1)-K*U(:,n-1));
      V(:,n-1)=(U(:,n)-U(:,n-2))/2/dt;
  end
  
  % figure,plot(t,F(1,:))
  figure,plot(t,U(1,:))
  figure,plot(t,U(2,:))
  figure,plot(t,F(1,:)+k*(U(2,:)-U(1,:)))
  
   %figure,plot([1:nstep-1]*dt,V(1,:))
   %figure,plot([1:nstep-1]*dt,AC(1,:))