Newmark法求振动方程，但是一直出现错误，这个程序是参考网上...

legendhu

clc;clear;close;
m=[106600,73500,7600,7600,73500,36600];
m=diag(m);    %质量矩阵
k=  [9.58*10^10 -3.15*10^10   0      0     0      0;
-3.15*10^10 9.43*10^10  -6.28*10^10   0     0      0;
0    -6.28*10^10 9.11*10^10  -2.83*10^10 0      0;
0     0     -2.83*10^10 9.11*10^10  -6.28*10^10  0;
0     0       0    -6.28*10^10 1.073*10^11  -4.45*10^10
0     0       0      0    -4.45*10^10 1.48*10^11];  %刚度矩阵
c=[4*exp(6) -4*exp(6) 0 0 0 0 ;
-4*exp(6) 4*exp(6) 0 0 0 0;
0 0 exp(6) -exp(6) 0 0;
0 0 -exp(6) exp(6) 0 0;
0 0 0 0 0 0;
0 0 0 0 0 0];    %阻尼矩阵
f0=1.589*10^7;
dt=0.0005;    %仿真步长
alfa=0.25;
beta=0.5;
tend=0.3
tt=0:dt:tend;

a0=1/alfa/dt/dt;
a1=beta/alfa/dt;
a2=1/alfa/dt;
a3=1/2/alfa-1;
a4=beta/alfa-1;
a5=dt/2*(beta/alfa-2);
a6=dt*(1-beta);
a7=dt*beta;
d=zeros(6,length(tt));
v=zeros(6,length(tt));
a=zeros(6,length(tt));
d(:,1)=[0,0,0,0,0,0]';
v(:,1)=[0,0,0,0,0,0]';
a(:,1)=[0,0,0,0,0,0]';
for i=1:length(tt);
f(:,i)=[f0*sin(500*tt),0,0,0,0,0]';
end

for i=2:length(tt);
ke=k+a0*m+a1*c;
fe(:,i-1)=f(:,i)+m*(a0*d(:,i-1)+a2*v(:,i-1)+a3*a(:,i-1))+c*(a1*d(:,i-1)+a4*v(:,i-1)+a5*a(:,i-1));
d(:,i)=inv(ke)*fe(:,i);
a(:,i)=a0*(d(:,i)-d(:,i-1))-a2*v(:,i-1)-a3*a(:,i-1);
v(:,i)=v(:,i-1)+a6*a(:,i-1)+a7*a(:,i);
end
plot(t,d)