how to do Newmark method for 2 different timesteps and to calculate mipoint deflection with respect to loads in all nodes?

Hi everyone,

i'm new to matlab.For my research, i am calculating dynamic analysis of beam for 2 different timesteps where i coulnt get the correct results. i have to calulate midpoint displacement of the beam when the load at each point (drawing attached). could anyone help me to solve this.

% Beam properties
L=10;                %length in m    length 10m 
A=0.00767;           %area in m^2
E=2.1e12;            %youngs modulus in N/m^2
ne1=10;              %no of elements
I=0.000636056;      %moment of inertia in m^4
nnp=ne1+1;           %no of node points
Le1=L\ne1;           %length of each element
M = 588.42;              % 60kg/m
P=10000;             %Force in 10kN 
ks = 4.5e+11;          %stiffness
gam=1/2;             
beta=1/4;            % average acceleration method  
C=0;                 %damping



% elementNodes: connections at elements
elementNodes=[1 2;2 3;3 4;4 5;5 6;6 7;7 8;8 9;9 10;10 11];
% numberElements: number of Elements
numberEements=size(elementNodes,1);
% numberNodes: number of nodes
numberNodes=11; 


%  elementNodes: connections at elements
 ii=1:numberElements;
 elementNodes(:,1)=ii;
 elementNodes(:,2)=ii+1;
 
%Dynamic load
% Time step

 t = linspace(0,11,11);  % 11 timestep and 110 timesteps (total simulation time 11sec)  
 dt = t(2)-t(1);        
 nt = length(t);


% Constants used in Newmark's integration
a1 = gam/(beta*dt) ; a2 = 1/(beta*dt^2) ;
a3 = 1/(beta*dt) ;       a4 = gam/beta ;
a5 = 1/(2*beta) ;        a6 = (gam/(2*beta)-1)*dt ;
depl = zeros(nt,1) ;
vel = zeros(nt,1) ;
accl = zeros(nt,1) ;
P = zeros(nt,1) ;

% Initial Conditions
depl(1) = 0 ;
depl(11) = 0;
vel(1) = 0 ;
P(1) = 10000;


accl(1) = (P(1) - C*vel(1) - ks*depl(1))/ M; 

kbar = ks + gam*C/(beta*dt) + M /(beta*dt*dt);
A = M /(beta*dt) + gam*C/beta;
B = M /(2*beta) +dt*C*((0.5*gam/beta)-1)*C;



% Time step starts

DPbar = zeros(nt,1);
Dv = zeros (nt,1);
veldv = zeros(nt,1);
accldv = zeros(nt,1);



%calculation of timestep
for i = 1:(length(t)-1)
     
DPbar(i) = P(1) + A*vel(i)+B*accl(i);
Dv(i) = DPbar(i)/kbar; 
depl(i+1) = depl(i) + Dv(i);
veldv(i) = gam*Dv(i)/(beta*dt) - gam*vel(i)/beta + dt*accl(i)*(1-0.5*gam/beta);
accldv(i)= Dv(i)/(beta*dt*dt) - vel(i)/(beta*dt) - accl(i)/(2*beta);
vel(i+1) = vel(i) + veldv(i);
accl(i+1) = accl(i) + accldv(i);
    
     
    
    
end