function [d,v,a]=Newmark(K,M,C,f,d1,v1,a1,dt,tend)
%K---剛度
%M---質量
%C---阻尼
%d1---初始位移
%v1---初始速度
%dt---時間步長
%tend---結束時間
%%%返回值
%d----位移
%v----速度
%a----加速度
gama=0.5;
beta=1/6;
alpha0=1/beta/dt^2;
alpha1=gama/beta/dt;
alpha2=1/beta/dt;
alpha3=1/2/beta-1;
alpha4=gama/beta-1;
alpha5=dt/2*(gama/beta-2);
alpha6=dt*(1-gama);
alpha7=gama*dt;
K1=K+alpha0*M+alpha1*C;
d=zeros(1,floor(tend/dt)+1);
v=zeros(1,floor(tend/dt)+1);
a=zeros(1,floor(tend/dt)+1);
d(1,1)=d1;
v(1,1)=v1;
a(1,1)=a1;%M\(f(1,1)-K*d1-C*v1);
for i=2:1:floor(tend/dt)+1
    t=(i-1)*dt;
    f1=f(1,i)+M*(alpha0*d(1,i-1)+alpha2*v(1,i-1)+alpha3*a(1,i-1))...
        +C*(alpha1*d(1,i-1)+alpha4*v(1,i-1)+alpha5*a(1,i-1));
    d(1,i)=K1\f1;
    a(:,i)=alpha0*(d(:,i)-d(:,i-1))-alpha2*v(:,i-1)-alpha3*a(:,i-1);
    v(:,i)=v(:,i-1)+alpha6*a(:,i-1)+alpha7*a(:,i);
end
return