[求助]请大家帮忙看看下列程序有错误吗

[chenjingqiang]

请大家帮忙看看下列程序有错误吗 运行时出现
??? Subscript indices must either be real positive integers or logicals.

请大家指教
%function [y2,y3,y4]=matr(Mass,Stiff,i,k)
Mass=[0.01 0.02 0.02;0.02 0.01 0.03;0.02 0.01 0.04];
Stiff=[2 3 4;5 4 6;2 1 3];
i=0.1;
k=0.2;
INVM=inv(Mass);
S=INVM*Stiff;
D=eig(S);
D1=sort(D); % 系统特征值从小到大
SQUAW1=D1(1); % 一阶特征值W1 的平方
fz=sqrt(SQUAW1)/2/3.14; % 一阶自振频率
A1=Stiff-SQUAW1*Mass;
V1=null(A1); % 对应一阶特征值的一阶正则振型向量
m01=V1'*Mass*V1;
k01=V1'*Stiff*V1;
ai=V1(i);% 系数
bk=V1(k);% 系数
A=[0 1;-SQUAW1 0];
B=[0;ai*bk/m01];
C=[1 0];D=[0];
Q=[0.001 0;0 0.001];% 给定控制误差的权值
R=[1]; % 给定控制能量的权值
K1=lqr(A,B,Q,R);
A2=A-B*K1;
%r=eig(A1);
[num1,den1]=ss2tf(A,B,C,D);
%[num2,den2]=ss2tf(A2,B,C,D);
% 保持系统在欠阻尼条件下 求取系统的超调量和调整时间 作为目标函数 外加系统一 阶自振频率
Wn=sqrt(den2(3)); % 系统自然振荡频率
b=den2(2)/Wn/2;
Qp=exp(-3.14*b/sqrt(1-b^2)); % 系统的超调量
Ts=(-log(0.05)-log(1-b^2)/2)/b/Wn; % 系统的调整时间衰减到0.02
% 建立多目标函数 一阶自振频率最大 超调量和调整时间最小
%f=3/fz+Ts+Qp/10
%输出数据
%y2=Qp;
%y3=Ts;
%y4=fz;
t=0:0.001:0.12;
%impulse(num2,den2,t); % 闭环系统冲激响应
step(num2,den2,t); % 闭环系统阶跃响应
%g=tf([num2],[den2]); % BODE 图
%bode(g,{10,100000});
%x0=[1 0.2];
%initial(A2,B,C,D,x0);% 连续系统的零输入响应

下面的是做了些改动的 是真正要用的请帮忙再看一下有问题吗 参数在前面已经生成,运行时y2,y3,y4返回值为零
function [y2,y3,y4]=matr(Mass,Stiff,i,k)
INVM=inv(Mass);
S=INVM*Stiff;
D=eig(S);
D1=sort(D); % 系统特征值从小到大
SQUAW1=D1(1); % 一阶特征值W1 的平方
fz=sqrt(SQUAW1)/2/3.14; % 一阶自振频率
A1=Stiff-SQUAW1*Mass;
V1=null(A1); % 对应一阶特征值的一阶正则振型向量
m01=V1'*Mass*V1;
k01=V1'*Stiff*V1;
ai=V1(i);% 系数
bk=V1(k);% 系数
A=[0 1;-SQUAW1 0];
B=[0;ai*bk/m01];
C=[1 0];D=[0];
Q=[0.001 0;0 0.001];% 给定控制误差的权值
R=[1]; % 给定控制能量的权值
K1=lqr(A,B,Q,R);
A2=A-B*K1;
%r=eig(A1);
[num1,den1]=ss2tf(A,B,C,D);
%[num2,den2]=ss2tf(A2,B,C,D);
% 保持系统在欠阻尼条件下 求取系统的超调量和调整时间 作为目标函数 外加系统一 阶自振频率
Wn=sqrt(den2(3)); % 系统自然振荡频率
b=den2(2)/Wn/2;
Qp=exp(-3.14*b/sqrt(1-b^2)); % 系统的超调量
Ts=(-log(0.05)-log(1-b^2)/2)/b/Wn; % 系统的调整时间衰减到0.02
% 建立多目标函数 一阶自振频率最大 超调量和调整时间最小
%f=3/fz+Ts+Qp/10
%输出数据
y2=Qp;
y3=Ts;
y4=fz;
t=0:0.001:0.12;
%impulse(num2,den2,t); % 闭环系统冲激响应
step(num2,den2,t); % 闭环系统阶跃响应
%g=tf([num2],[den2]); % BODE 图
%bode(g,{10,100000});
%x0=[1 0.2];
%initial(A2,B,C,D,x0);% 连续系统的零输入响应