how to perform Newmark method for 2 different time step and to find defection at one point with respect to load on all nodes?

Hi everybody,

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 is 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
L e1 = L \ n e1;           % length of each element
M = 588.42;              % 60kg / m
P = 10,000;             % 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
numberElments = size (element nodes, 1);
% numberNodes: number of nodes
number nodes = 11; 


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

t = linspace (0,11,10);  % 11 timesteps and 110 timesteps (total simulation time 11sec)  
 dt = t (2) -t (1);
n = length (m);         % no of nodes
 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