function [x,xdot,x2dot]=newmark_B(K,C,M,beta,gamma,x0,xdot0,x2dot0,delt,t,ft) a0=1/(beta*delt^2) a1=gamma/(beta*delt) a2=1/(beta*delt) a3=(.5/beta)-1 a4=(gamma/beta)-1 a5=.5*delt*((gamma/beta)-2) a6=delt*(1-gamma) a7=gamma*delt Keff=K+a0*M+a1*C A=choleski(Keff) l=0:delt:t; k1=length(l) x(:,1)=x0;xdot(:,1)=xdot0;x2dot(:,1)=x2dot0; for i=2:k1 feff=ft(:,i)+M*(a0*x(:,(i-1))+a2*xdot(:,(i-1))+a3*x2dot(:,(i-1)))+C*(a1*x(:,(i-1))+a4*xdot(:,(i-1))+a5*x2dot(:,(i-1))); x(:,i)=subsitution(A,feff); x2dot(:,i)=a0*(x(:,(i))-x(:,(i-1)))-a2*xdot(:,(i-1))-a3*x2dot(:,(i-1)); xdot(:,i)=xdot(:,(i-1))+a6*x2dot(:,(i-1))+a7*x2dot(:,(i)); end