This is for single degree of freedom but when I am writing it for MDOF why the DFbar is not calculating properly?

function [t,u,v,a, kbar, R, keyp, key, Du, k_p] = NewmarkNon( t, DF, Dt, u0, v0, ut, uc,a0, F1, k, c, m, Rt, Rc, R, gamma, beta, key, k_p)

kbar = k_p +gamma*c/(beta*Dt)+m/(beta*Dt*Dt);
A = m/(beta*Dt)+gamma*c/beta; 
B = m/(2*beta)+Dt*c*((0.5*gamma/beta)-1); 
DFbar = DF + A*v0+B*a0; % Incremental effective force (Eq.6.59)    I am not getting DFbar=DF for i=1, when v0=0 and a0=0 for i=1
Du = DFbar/kbar; 
Dudot = gamma*Du/(beta*Dt)-gamma*v0/beta+ Dt*a0*(1-0.5*gamma/beta); 
u=u0+Du ; 
v=v0+Dudot; 
a=1/m*(F1-c*v-R);
t=t+Dt;
end