how to perform Newmark method for fixed beam?

hi,

i am new to matlab. i am doing dynamic analysis of a beam using newmark method with different time steps at the midpoint point of node 6 of the beam. the outputs are not correct as i expected. Displacement is linear and not parabloic. could anyone help me out where did i do the mistake.

% 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
n=10;              %no of elements
I=0.000636056;      %moment of inertia in m^4
nn=n+1;           %no of node points
Le=L\n;           %length of each element
m = 588.42;              % 60kg/m
P=10000;             %Force in 10kN 
p = 7800;
%ks = 4.5e+11;          %stiffness
Gamma=1/2;             
Beta=1/4;            % average acceleration method  
C=0;                 %damping


%global stiffness  and mass matrix for fixed beam

Ka =  E*I/Le^3 *[24 -12 0 0 0 0 0 0 0;
                 -12 24 -12 0 0 0 0 0 0;
                 0 -12 12 -12 0 0 0 0 0;
                 0 0 -12 12 -12 0 0 0 0;
                 0 0 0 -12 12 -12 0 0 0;
                 0 0 0 0 -12 12 -12 0 0;
                 0 0 0 0 0 -12 12 -12 0;
                 0 0 0 0 0 0 -12 24 -12;
                 0 0 0 0 0 0 0 -12 24];


Ma =  p*A*Le/2 * [588.42 0 0 0 0 0 0 0 0;
                  0 588.42 0 0 0 0 0 0 0;
                  0 0 588.42 0 0 0 0 0 0;
                  0 0 0 588.42 0 0 0 0 0;
                  0 0 0 0 588.42 0 0 0 0;
                  0 0 0 0 0 588.42 0 0 0;
                  0 0 0 0 0 0 588.42 0 0;
                  0 0 0 0 0 0 0 588.42 0;
                  0 0 0 0 0 0 0 0 588.42];

%C = Dalpha * Ka + Dbeta * Ma;

% integration constant
c1 = 1/Beta/dt^2;
c2 = Gamma/Beta/dt;
c3 = 1/Beta/dt;
c4 = 1/2/Beta-1;
c5 = Gamma/Beta-1;
c6 = (Gamma/2/Beta-1)*dt;
c7 = (1-Gamma)*dt;
c8 = Gamma*dt;

% % Time step

 t = linspace(1,10,11);  % 10 timestep 
 dt = 1;
 nt = length(t);
 
 
 P =zeros(9);
 P(:,1) = 10000;

 
 %intial conditions

Fbar = zeros(9);
Dv = zeros(9);
veldv = zeros(9);
accldv = zeros(9);
D = zeros(9);
V = zeros(9);
Aa = zeros(9);
 D(:,1) = 0;
 V(:,1) = 0;
 Aa(:,1) = (Ma^-1)*(P(:,1) - C*V(:,1) - Ka*D(:,1));
 
Kbar = c1*Ma+c2*C+Ka;
%A = Ma /(Beta*dt) + Gamma*C/beta;
%B = Ma /(2*Beta) +dt*C*((0.5*Gamma/Beta)-1)*C;

for i = 1:(nt-1)

    Fbar(:,i+1) = P(:,1) + Ma * (c1 * D(:,i) + c3 * V(:,i) + c4 * Aa(:,i)) + C * (c2 * D(:,i) + c5 * V(:,i) + c6 * Aa(:,i));  
    Dv(:,i+1) = Kbar^-1 * Fbar(:,i+1);
    D(:,i+1) = D(:,i) + Dv(:,i+1);    
    accldv(:,i+1) = ((c1 * D(:,i+1)) - D(:,i)) - (c3 *V(:,i)) - (c4 * Aa(:,i));
    Aa(:,i+1) = Aa(:,i) + accldv(:,i+1);
    veldv(:,i+1) = V(:,i) + (c7 * Aa(:,i)) + (c8 * Aa(:,i+1));
    V(:,i+1) = V(:,i) + veldv(:,i+1);   
    
end