% CENTRAL DIFFERENCE
% femem_springdyn_cd1
% RAMP 1DOF
clear
   nstep=1000;
   T=4;
   T0=0.2;
   F0=2000;
   M=31.83;
   K=100;
    
   dt=T/nstep;
   t=[1:nstep]*dt;
   t1=[1:T0/dt]*dt;
   t2=[T0/dt+1:T/dt]*dt;
   
   NumEQ=1;
   
   F=zeros(NumEQ,nstep);
   V=zeros(NumEQ,nstep);
   U=zeros(NumEQ,nstep);
   AC=zeros(NumEQ,nstep);
   
   F(1:T0/dt+1)=F0*(1-dt/T0*[0:T0/dt]);
  
   
   AC(:,1)=(F(:,1)-K*U(:,1))/M;
   UP=U(:,1)-dt*V(:,1)+dt^2/2*AC(:,1);
   U(:,2)=(dt^2*F(:,1)+(2*M-dt^2*K)*U(:,1)-M*UP)/M;
   
  for n=3:nstep
      U(:,n)=(dt^2*F(:,n-1)+(2*M-dt^2*K)*U(:,n-1)-M*U(:,n-2))/M;
      AC(:,n-1)=(F(:,n-1)-K*U(:,n-1))/M;
      V(:,n-1)=(U(:,n)-U(:,n-2))/2/dt;
  end
  
  w=(K/M)^.5;
  
  Uexact(1:T0/dt)=F0/K*(1-cos(w*t1))+F0/K/T0*(sin(w*t1)/w-t1);
  Uexact(1+T0/dt:T/dt)=F0/K/T0/w*(1-cos(w*T0))*sin(w*t2)+F0/K*(sin(w*T0)/w/T0-1)*cos(w*t2);
  
  Vexact(1:T0/dt)=F0/K*w*sin(w*t1)+F0/K/T0*(cos(w*t1)-1);
  Vexact(1+T0/dt:T/dt)=F0/K/T0*(1-cos(w*T0))*cos(w*[1+T0/dt:T/dt]*dt)-F0/K*(sin(w*T0)/w/T0-1)*w*sin(w*[1+T0/dt:T/dt]*dt);
  
  ACexact(1:T0/dt)=F0/K*w*w*cos(w*t1)-w*F0/K/T0*sin(w*t1);
  ACexact(1+T0/dt:T/dt)=-w*F0/K/T0*(1-cos(w*T0))*sin(w*[1+T0/dt:T/dt]*dt)-F0/K*(sin(w*T0)/w/T0-1)*w*w*cos(w*[1+T0/dt:T/dt]*dt);
  
  figure,plot([1:nstep],U,'-',[1:nstep],Uexact)
  figure,plot([1:nstep],V,'-',[1:nstep],Vexact)
  figure,plot([1:nstep],AC,'-',[1:nstep],ACexact)