求救，matlab单自由度的时程分析程序

ultiwin

如题，自己按照Newmark-beta法编制了一个单自由度的时程分析程序，但算出来的结果和sap2000的不符，而且反应加速度值很大，不知道问题出在哪里，源程序如下：

clear;clc;
% 结构模型初始参数----------------------------------------------------------
m=3e3;                    %质量(单位：kg)
k=1e6;                    %刚度（(单位：N/m)）
kesai=0.05;               %阻尼比取0.05
c=2*kesai*sqrt(k*m);      %阻尼系数
% 读取地震波数据------------------------------------------------------------
acc=textread('D:\研究生课题资料\科研课题资料\处理后的smc文件\51WCW_90_chnua370295.smc_090501.a','%f','headerlines',56);
PGA_Max=max(abs(acc))     %最大地面加速度绝对值
% Newmark-beta法的基本参数--------------------------------------------------
beta=1/6;  gama=0.5;      %按线性加速度法计算更接近真实结果，故取此组参数
dt=0.02;                  %地震加速度时程波记录时间间隔
b1=1/(beta*dt^2);     b2=1/(beta*dt);     b3=1-1/(2*beta);    %计算参数
b4=gama/(beta*dt);    b5=gama/beta-1;     b6=(1-gama/(2*beta))*dt;
ke=k+m*b1+c*b4;           %等效刚度
% 设定结构初始状态为零,生成向量空间存储计算值---------------------------------
u=zeros(100/dt,1);   v=zeros(100/dt,1);   a=zeros(100/dt,1);  
% Newmark-beta法的主计算程序------------------------------------------------
for n=2:100/dt
    fe=-m*acc(n)+[b1*u(n-1)+b2*v(n-1)-b3*a(n-1)]*m+[b4*u(n-1)+b5*v(n-1)-b6*a(n-1)]*c;  %等效荷载   
    u(n)=fe/ke;
    a(n)=b1*[u(n)-u(n-1)]-b2*v(n-1)+b3*a(n-1);
    v(n)=b4*[u(n)-u(n-1)]-b5*v(n-1)+b6*a(n-1);  
end
% 绘制结构在地震作用下的位移、速度、加速度时程曲线-----------------------------
subplot(3,1,1)
t=(0:length(a)-1)*dt;
plot(t,a)    %加速度时程曲线
Acc_Max=max(abs(a))
title('Earthquake Response Curve Of Station 51WCW-90','fontsize',15)
ylabel('Acc(cm/s^2)','fontsize',12)
subplot(3,1,2)
plot(t,v)    %速度时程曲线
Vel_Max=max(abs(v))
ylabel('Vel(cm/s)','fontsize',12)
subplot(3,1,3)
plot(t,u)    %位移时程曲线
Dis_Max=max(abs(u))
xlabel('Time/s','fontsize',12)
ylabel('Dis/cm','fontsize',12)
% End---程序结束-------------
小弟初次发贴，恳请达人们帮分析一下，不胜感激！

ultiwin

其中的循环部分是根据结构动力学书上的写的，感觉问题就出在那部分了，请高人们指点一下

ultiwin

急啊，请达人们帮下忙

ultiwin

附地震波数据文件

ultiwin

原文件.a格式的，传不上来，所以改成txt格式了。期待达人们的解答