clear
clc

A=1 ;
nn=10;
M=zeros(nn,1);%global mass matrix
K=zeros(nn,nn);

Mel= [10,0;0,10];% element mass
Kel= [1100,-1100;-1100,1100]; % element stiffness matrix

for I=2:nn-1
    M(I,I)=10;
end
Kel=[1100,-1100;-1100,1100];
for i=1:nn-1
     j=i+1;

K=LinearBarAssemble (K,Kel,i,j);
end

%time step
dt_crit=2*sqrt (10/1100);

factor=.009;
dt=dt_crit*factor;
%dt=0.4;

%total time
t=2;
tn=round (t/dt);% number of total time step

%initial conditions
d0= 0;
v0= A*sqrt (1100/10); %
a0= 0;

%define vector
a=zeros (nn,1);
v=zeros (nn,1);
d=zeros (nn,1);

%load
%Rext=zeros (nn,1);
Rint=zeros (nn,1);
Rext=zeros(tn,1);

M(1)=5;
M(2:nn-1)=10;
M(nn)=5;

%start
d(1)= 0;
v(1)= 0;
d(2)= d0+dt*v0+(0.5*dt^2)/M(1)*(Rext(1)-Rint(1));  % lecture4&5 slide 24&23
%d(2)= d0+dt*v0+(0.5*dt^2)*M^(-1)*(100-0);  % lecture4&5 slide 24&23

v(2)= v0+0.5*dt*a0;

v_time=zeros(tn,1);
v_d=zeros(nn,tn);

v_node=[1,2,3,4,5,6,7,8,9,10];

for i=2:tn-1
  
     if i<300
        Rext(i)=100;
     else
        Rext(i)=0;
     end
  for ni=1:nn
    a(ni)=(1/M(ni))*(Rext(ni)-Rint(ni));
    d(ni+1)=d(ni)+dt*v(ni)+dt^2*a(ni);
    v(ni+1)=v(ni)+dt*a(ni);
    v_time(i,1)=i*dt;
    
    v_d (ni,i)=d(ni);
  end
    
  for ns=1:nn-1%stress 
      f=1*200*(d(ni+1)-d(ni))/2;% f=A*e*strain
      Rint(ns)=f+Rint(ns);
      Rint(ns+1)=Rint(ns+1)-f;
  end 

    %v_d(1:nn,i-1)=d(ni);
    mesh(v_time,v_node,v_d)
end


%v_node=[1,2,3,4,5,6,7,8,9,10];
%mesh(v_time,v_node,v_d)