请问各位大牛怎么用matlab编程求系统的频谱特性曲线和时程曲线

ME!

function [d,v,a] = Wilson( K, M, C, f, d1, v1, dt, tend )
% 利用Wilson-theta 法计算结构的动力响应
K=sysK;% K ----- 调用刚度矩阵
M=sysM; % M ----- 调用质量矩阵
  C=0;
dt=0.1;% dt ----- 时间步长
tend=40;% tend --- 结束时间
w=600;
[n,n] = size( K ) ;

%Wilson算法
theta = 1.37 ;
tao = theta*dt ;
alpha0 = 6/tao^2 ;
alpha1 = 3/tao ;
alpha2 = 2*alpha1 ;
alpha3 = tao/2 ;
alpha4 = alpha0/theta ;
alpha5 = -alpha2/theta ;
alpha6 = 1-3/theta ;
alpha7 = dt/2 ;
alpha8 = dt^2/6 ;
K1 = K + alpha0*M + alpha1*C ;
d = zeros( n, floor(tend/dt) + 1 ) ;
v = zeros( n, floor(tend/dt) + 1 ) ;
a = zeros( n, floor(tend/dt) + 1 ) ;
f = zeros( n, floor(tend/dt) + 1 ) ;
d(:,1) = 0.01 ; %初始位移
  v(:,1) = 0 ;%初始速度
f(:,1) =0 ;%初始载荷
a(:,1) = inv(M)*(f(:,1)-K*d(:,1)-C*v(:,1)) ;  %初始加速度
t=0:dt:tend;
for i=2:1:length(t)
%t(i) = (i-1)*dt ;
f(:,i)=0*100*sin(w*t(i));
ftheta = floor(theta) ;
%fq = f(i-1+ftheta-1)+ (theta-ftheta)*( f(i+ftheta-1) - f(i+ftheta-2) ) ;
fq=f(:,i-1)+ftheta*(f(:,i)-f(:,i-1));
f1 = fq + M*(alpha0*d(:,i-1)+alpha2*v(:,i-1)+2*a(:,i-1))+ C*(alpha1*d(:,i-1)+2*v(:,i-1)+alpha3*a(:,i-1)) ;
dq= inv(K1)*f1 ;
a(:,i) = alpha4*(dq-d(:,i-1)) + alpha5*v(:,i-1) + alpha6*a(:,i-1) ;
v(:,i) = v(:,i-1) + alpha7 * ( a(:,i) + a(:,i-1) ) ;
d(:,i) = d(:,i-1) + dt*v(:,i-1) + alpha8 * ( a(:,i)+2*a(:,i-1) ) ;
end

  
    %Newmark算法
gama = 0.5 ;
beta = 0.25 ;
[n,n] = size( K ) ;
Nalpha0 = 1/beta/dt^2 ;
Nalpha1 = gama/beta/dt ;
Nalpha2 = 1/beta/dt ;
Nalpha3 = 1/2/beta - 1 ;
Nalpha4 = gama/beta - 1 ;
Nalpha5 = dt/2*(gama/beta-2) ;
Nalpha6 = dt*(1-gama) ;
Nalpha7 = gama*dt ;
NK1 = K + Nalpha0*M + Nalpha1*C ;
Nd = zeros( n, floor(tend/dt) + 1 ) ;
Nv = zeros( n, floor(tend/dt) + 1 ) ;
Na = zeros( n, floor(tend/dt) + 1 ) ;
Nd(:,1) = 0.01; %初始位移
  f(:,1) =0 ;%初始载荷
Na(:,1) = inv(M)*(f(:,1)-K*Nd(:,1)-C*Nv(:,1)) ;  %初始加速度
t=0:dt:tend;
for i=2:1:length(t)
f(:,i)=0*100*sin(w*t(i));
f2 = f(:,i) + M*(Nalpha0*Nd(:,i-1)+Nalpha2*Nv(:,i-1)+Nalpha3*Na(:,i-1))+ C*(Nalpha1*Nd(:,i-1)+Nalpha4*Nv(:,i-1)+Nalpha5*Na(:,i-1)) ;
Nd(:,i) = inv(NK1)*f2 ;
Na(:,i) = Nalpha0*(Nd(:,i)-Nd(:,i-1)) - Nalpha2*Nv(:,i-1) - Nalpha3*Na(:,i-1) ;
Nv(:,i) = Nv(:,i-1) + Nalpha6*Na(:,i-1) + Nalpha7*Na(:,i) ;
end
plot(t,d(1,:),'-b^',t,Nd(1,:),'g',t,x,'r')
xlabel('t');
ylabel('u(t)');
title('位移与时间的关系');
还有就是那个载荷，我是给的0，不知道怎么加
我的模型是实体的，总体质量矩阵和刚度矩阵均为51阶的方阵
我验证了wilson和Newmark方法对单自由度和两自由度的时间和位移的响应曲线与解析解的结果一致
但是对于我的模型，我画出来的曲线却是一条近似为斜率为1的直线，我不懂是为什么