程序在下面，运行时提示“Matrix mast be square”，按照网上有人说的把“*”全部改为“.*”,还是报同样的错误
[code]disp('**************************************************').
disp('*** 基于抛物线内插法的不同阻尼比情况下反应谱计算  ***')
disp('***          输入地震加速度为X``(t)=cos2πt       ***')
disp('***         Author:WanLi Zhao  (BENG)          ***')
disp('***           Email:ewarren@126.com            ***')
disp('***             HOHAI  University              ***')
disp('***           2010.05.13  NanJing              ***')
disp('**************************************************')
g=('请输入阻尼比：');

  for i=1:1:1000
    for k=0:1:1000
        j=2*pi;            %θ=2π
        w=(2*pi)/(i/100);    %自振圆频率ω
        t=k/100;           %时间
        u(1)=0;            %初始位移为0
        u1(1)=0;           %初始速度为0
        c11=exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01+(g/sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01;
        c12=exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01;
        d11=(1/w^2)*(1+3*g/w*0.01+(4*g^2-1)/w^2*0.01^2)*exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01+(g+3*g/w*0.01+g(4*g^2-1)/w^2*0.01^2-3/2*w*0.01-2*g/w^2*0.01^2)*(1/w^2*sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01-(1/w^2)*(g/w*0.01+(4*g^2-1)/w^2*0.01^2);
        d12=-(2/w^3*0.01)*(2*g+(4*g^2-1)/w*0.01)*exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01-(2*g^2+(4*g^2-1)*g/w*0.01-1-2*g/w*0.01)*(2/0.01*w^3*sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01+(1/w^2*0.01)*(2*(4*g^2-1)/w^2*0.01-0.01);
        d13=(1/w^3*0.01)*(g+(4*g^2-1)/w^2*0.01)*exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01+(g^2+(4*g^2-1)*g/w*0.01-1/2-2*g/w*0.01)*(1/0.01*w^3*sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01-((4*g^2-1)/w^2*0.01-g/w)*(1/w^2*0.01);
        c21=-(w/sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01;
        c22=-(g/sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01+exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01;
        d21=-(3+4*g/w*0.01)*(1/2*w^2*0.01)*exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01-(1+3*g/2*w*0.01+(2*g^2-1)/w^2*0.01^2)*(1/w*sqrt(1-g^2));
        d22=(1+2*g/w*0.01)*(2/w^2*0.01)*exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01+(g+(2*g^2-1)/w*0.01)*(2/0.01*w^2*sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01-4*g/w^3*0.01^2;
        d23=-(1+4*g/w*0.01)*(1/2*w^2*0.01)*exp(-g*w*0.01)*cos(w)*sqrt(1-g^2)*0.01-(g/2+(2*g^2-1)/w*0.01)*(1/0.01*w^2*sqrt(1-g^2))*exp(-g*w*0.01)*sin(w)*sqrt(1-g^2)*0.01+(1/(2*w^2*0.01))*(4*g/w*0.01-1);
        u(k+2)=c11*u(k)+c12*u1(k)+d11*cos(j*(k+1))+d12*cos(j*(k+2))+d13*cos(j*(k+3));
        u1(k+2)=c21*u(k)+c22*u1(k)+d21*cos(j*(k+1))+d22*cos(j*(k+2))+d23*cos(j*(k+3));
        u11(k+2)=-2*g*w*u1(k+1)-w^2*u(k+1);
    end
    wy=max(u);
    sd=max(u1);
    jsd=max(u11);
end    
         %-------输出结果-------
    fid=fopen('抛物线内插法位移解.txt','w');
    fprintf(fid,'%15.5f\r\n',wy);
    fclose(fid);
    fid=fopen('抛物线内插法速度解.txt','w');
    fprintf(fid,'%15.5f\r\n',sd);
    fclose(fid);
    fid=fopen('抛物线内插法加速度解.txt','w');
    fprintf(fid,'%15.5f\r\n',jsd);
    fclose(fid);

    %---------画图-------
    figure;
    T=0.01:0.01:10;
    plot(T,wy);
    title('y(t)');
    xlabel('T');
    ylabel(g);
    grid on;
    %--------画图-------
    figure;
    T=0.01:0.01:10;
    plot(T,sd);
    title('y′(t)');
    xlabel('T');
    ylabel(g);
    grid on;
    %--------画图--------
    figure;
    T=0.01:0.01:10;
    plot(T,jsd);
    title('y′′(t)');
    xlabel('T');
    ylabel(g);
    grid on;
        [code]