使用传递函数求解响应的程序问题！！

SaintYZ

使用下面的程序可以算出结果，可是结果却不对，频响函数法我想各位大牛肯定很熟了，不知程序那部分写法有错，还是理论上有问题？？？
望指点！！！
××××××程序是3个自由度的小车系统，每个小车上作用一个与质量成比例的力，求响应×××××××
function  y=transfer   
m1=[1 1 1];
m=diag(m1);                %形成3×3的质量阵
k=[3 -1 0;-1 2 -1;0 -1 3]; %形成3×3的刚度阵
[a,b]=eig(k,m);
d=diag(sqrt(b));           %求解固有频率和固有振型
for i=1:3
[d1(i),j]=min(d);
xgd(:,i)=a(:,j);
d(j)=max(d)+1;
end
w=d1;
u=xgd;                       %将固有频率和相应的振型按升序排列
v=u';                       %解耦质量阵，刚度阵
m0=v*m*u;
k0=v*k*u;
syms w t;
for p=1:3                    %求频响传递函数矩阵
    for k=1:3
        h1(p,k)=u(p,1)*u(k,1)/(k0(1,1)-w.^2*m0(1,1));
        h2(p,k)=u(p,2)*u(k,2)/(k0(2,2)-w.^2*m0(2,2));
        h3(p,k)=u(p,3)*u(k,3)/(k0(3,3)-w.^2*m0(3,3));
        H(p,k)=h1(p,k)+h2(p,k)+h3(p,k);
    end
end
HH=-w.^2*H;    %加速度传递函数矩阵
f=sin(6*pi*t)*m1;    %力激励

for i=1:3
    F(i)=fourier(f(i),t,w);     %把激励转换成频率域
  
end
X=F*H;
Y=F*HH;                         %计算频域中的位移和加速度
for i=1:3
    x(i)=ifourier(X(i),w,t);
    y(i)=ifourier(Y(i),w,t);
end                          %计算时域中的位移和加速度
t=0.006;                     %0.025时刻的加速度数值解
xx1=real(subs(x));
yy1=real(subs(y));
disp(xx1);
disp(yy1);