matlab请问这个程序哪里错了

matlab请问这个程序哪里错了
%主程序（matlab）：
x=[0.137 65*pi/180 0.1 0.125 0.6 0.15 0.275 0.575 0];%赋初值
dr=pi/180;%度转化为弧度
dth=10*dr;w1=1;%每10度计算一个点
for i=1:37
y=ntpc(x);%调用从动件位置方程求解函数ntpc(自编)
s3=y(1);
theta3=y(2);
theat4=y(4);
se=y(4);%得到位置参数.
%将各位置参数用向量储存,便于后面绘图,角度用度表示
ss3(i)=y(1);th1(i)=x(10)/dr;th3(i)=y(2)/dr;th4(i)=y(3)/dr;sse(i)=y(4);
%进行速度分析
A=[cos(theta3) -s3*sin(theta3) 0 0;
sin(theta3) s3*cos(theta3) 0 0;
0 -x(6)*sin(theta3) -x(7)*sin(theta4) -1;
0 x(6)*cos(theta3) x(7)*cos(theta4) 0];%A机构从动件的位置参数矩阵
B=[-x(5)*sin(x(10));x(5)*cos(x(10));0;0];%B机构原动件的位置参数列阵
yy=w1*inv(A)*B;%公式1-3求解,yy表示机构从动件速度列阵,inv(A)是A的逆阵
vs3=yy(1);w3=yy(2);w4=yy(3);vse=yy(4);
%将各速度参数以向量的方式表示,以便后面绘图
dvs3(i)=yy(1);dw3(i)=yy(2);dw4(i)=yy(3);dvse(i)=yy(4);
%dA为从动件位置参数矩阵对时间一次求导
%进行角速度分析
dA=[-w3*sin(theta3) -vs3*sin(theta3)-s3*w3*cos(theta3) 0 0;
w3*cos(theta3) vs3*cos(theta3)-s3*w3*sin(theta3) 0 0;
0 -x(6)*w3*cos*(theta3) -x(7)*w4*cos(theta4) 0;
0 -x(6)*w3*sin*(theta3) -x(7)*w4*sin(theta4) 0];
%dB就是原动件位置参数列阵对时间一次求导
dB=[-x(5)*w1*cos(x(10));-x(5)*w1*sin(x(10));0;0];
KK=-dA*yy+w1*dB;%KK为公式1-4右端
ya=inv(A)*KK;%公式1-4求解,ya为从动件加速度列阵
%将各加速度以向量表示
as3(i)=ya(1);atheta3(i)=ya(2);atheta4(i)=ya(3);ase(i)=ya(4);
x(10)=x(10)+dth;%计算下一个点
x(1)=s3;
x(2)=theta3;
x(3)=theta4;
x(4)=se;
end
%绘制运动参数曲线
ssubplot(2,2,1); %选择第1个子窗口
plot(th1,th3,th1,th4,th1,sse*1e3);%绘制位置线图
subplot(2,2,2);
plot(th1,dw3,th1,dw4,th1,dvse);%绘制速度线图
subplot(2,2,3);
plot(th1,atheta3,th1,atheta4,th1,ase);%绘制加速度线图
%这个函数是关于牛头刨床位置方程求解,可得：s3,theta3,theta4,sE四个运动变量.
%子程序：
function[y]=ntpc(x)
s3=x(1);
theta3=x(2);
theta4=x(3);
se=x(4);
epsilon=1e-6;%设置求解精度为10-6
%用矩阵的形式表示位置方程组（4x1的矩阵）
f=[s3*cos(theta3)-x(5)*cos(x(10));
s3*sin(theta3)-x(5)*sin(x(10))-x(8);
x(6)*cos(theta3)+x(7)*cos(theta4)-se;
x(6)*sin(theta3)+x(7)*sin(theta4)-x(9)];
%用牛顿-辛普森法求解
while norm(f)>epsilon
%J位置方程组的雅可比矩阵,即从动件位置参数矩阵
J=[cos(theta3) -s3*sin(theta3) 0 0;
sin(theta3) s3*cos(theta3) 0 0;
0 -x(6)*sin(theta3) -x(7)*sin(theta4) -1;
0 x(6)*cos(theta3) x(7)*sin(theta4) 0];
dth=inv(J)*(-1.0*f);%计算增量,进行迭代,inv(J) 为J的逆阵
s3=s3+dth(1);
theta3=theta3+dth(2);
theta4=theta4+dth(3);
se=se+dth(4);
f=[s3*cos(theta3)-x(5)*cos(x(10));s3*sin(theta3)-x(5)*sin(x(10))-x(8);
x(6)*cos(theta3)+x(7)*cos(theta4)-se;x(6)*sin(theta3)+x(7)*sin(theta4)-x(9)];
norm(f);%若未达精度,会继续迭代.
end
%输出4个参数
y(1)=s3;
y(2)=theta3;
y(3)=theta4;
y(4)=se;

我不知道你的算法所以没法帮你改,但是我可以从语法上帮你指出哪里错误了
所先从%子程序：
function[y]=ntpc(x)
到结尾y(4)=se;这些代码你要剪切下来另存为ntpc.m
其次
y=ntpc(x);%调用从动件位置方程求解函数ntpc(自编)
在调用ntpc(x)时
f=[s3*cos(theta3)-x(5)*cos(x(10));
x(10)就溢出了.
自己看看吧
再问： 溢出了该怎么改，我不太会
再答： x=[0.137 65*pi/180 0.1 0.125 0.6 0.15 0.275 0.575 0] 有9个数据，再添加一个吧，看看是不是少了 我也不知 你的程序算法原理，无从改起啊。