function [T,q]=Newmark3(M,B,C,K,f,q0,qd0,qd20,alpha,beta,gamma,h,Tf)
syms t
N=size(M,1);
T=0:h:Tf;
lt=length(T);
F=subs(f,t,T);
qd30=M\(F(1)-B*qd20-C*qd0-K*q0);
%===
q=zeros(N,lt);
qd=zeros(N,lt);
qd2=zeros(N,lt);
qd3=zeros(N,lt);
q(:,1)=q0;
qd(:,1)=qd0;
qd2(:,1)=qd20;
qd3(:,1)=qd30;
%===
A=M+alpha*h*B+gamma*(h^2)*C+beta*(h^3)*K;
for n=1:lt-1
    D=F(n+1)-B*(qd2(:,n)+(1-alpha)*h*qd3(:,n))-...
        C*(qd(:,n)+h*qd2(:,n)+(0.5-gamma)*(h^2)*qd3(:,n))-...
        K*(q(:,n)+h*qd(:,n)+(h^2)/2*qd2(:,n)+(1/6-beta)*(h^3)*qd3(:,n));
    qd3(:,n+1)=A\D;
    qd2(:,n+1)=qd2(:,n)+(1-alpha)*h*qd3(:,n)+alpha*h*qd3(:,n+1);
    qd(:,n+1)=qd(:,n)+h*qd2(:,n)+(0.5-gamma)*(h^2)*qd3(:,n)+gamma*(h^2)*qd3(:,n+1);
    q(:,n+1)=q(:,n)+h*qd(:,n)+(h^2)/2*qd2(:,n)+(1/6-beta)*(h^3)*qd3(:,n)+...
        beta*(h^3)*qd3(:,n+1);
end