请问问题在哪，检查了好几遍，求解释
main.m
%问题1
Lc = 22.02;%锚链的长度
cId = 2;%锚链的型号
M = 1200;%重物球的质量
depth = 18;%海的深度
w = 12;%w风的速度
v = 0;%v水的速度
%tilt:钢桶与竖直方向的角度
%elev:最后一根杆与海底平面的角度(锚链跟海底平面的夹角)
%xsbed:躺在海底的锚链，最大的是躺在海底的最后一节链节
%xbuoy:最后一根杆的x值，即与浮标相连的钢管的x的值就是浮标的游动距离
%f:吃水深度与整个浮标的比例
[tilt,elev,xsbed,xbuoy,f] = moor(Lc,cId,w,v,M,depth)
moor.h
%Lc锚链的长度，cId锚链的型号，w风的速度，v水的速度，M重物球的质量，depth海的深度
function [tilt,elev,xsbed,xbuoy,f] = moor(Lc,cId,w,v,M,depth)
%g重力加速度 rho海水的密度，rhoFe铁的密度,cdwin风的阻力系数
g = 9.81;rho = 1.025e3;rhoFe = 7.9e3;cdwin = 0.625;
%lc一节锚链的长度，mc一节锚链的质量，dc一节锚链的直径，Lc锚链总长度,nc锚链的节数
[lc,mc,dc] = chainpara(cId);nc = round(Lc/lc);
%m所有杆的质量，h所有杆的长度，d杆的直径
m = [1000,10*ones(1,4),100,mc*ones(1,nc)];
h = [2, ones(1,4), 1, lc*ones(1,nc)];
d = [2, 5e-2*ones(1,4),0.3,dc*ones(1,nc)];
Fb = pi*(d/2).^2.*h*rho*g - m*g;%Fb所有杆的浮力，浮标除外
phi = zeros(1,length(h));%每一根杆的倾斜角度初始化为0
fmin = 0;fmax = 1;%设置吃水比例，最小值和最大值
while fmax - fmin > 1e-10
f = (fmax + fmin)/2;
Fb(1) = rho*f*pi*(d(1)/2).^2.*h(1)*g - m(1)*g;%浮标的浮力
Fw = cdwin * (1-f)*h(1).*d(1)* w .^2;%风的阻力
Fs = waterload(v,h,d,phi,depth,f);%水的作用力，水上(水下可以不考虑)
phi = solvequileq(Fb,Fw,Fs,M,f);%计算出杆的倾斜角度
x = h.*sind(phi);%x杆在水平方向的投影
y = h.*cosd(phi);%y杆在竖直方向上的投影
if sum(y(2:end)) + h(1)*f>depth; fmax = f;else;fmin = f;end
end
x = cumsum([0 flip(x)]);y = cumsum([0 flip(y)]);%所有杆的x和y
tilt = phi(6); %钢桶与竖直方向的角度
elev = 90 - phi(end);%最后一根杆与海底平面的角度(锚链跟海底平面的夹角)
xsbed = max(x(y<1e-10));%躺在海底的锚链，最大的是躺在海底的最后一节链节
xbuoy = x(end - 1);%最后一根杆的x值，即与浮标相连的钢管的x的值就是浮标的游动距离
waterload.m
%v:水的速度,h：杆的长度,d：杆的直径,phi：杆与竖直方向的角度,depth：水的深度
%f：吃水深度与总高度的比例
%Fs水的作用力，水上(水下可以不考虑)
function Fs = waterload(v,h,d,phi,depth,f)
y = h.*cosd(phi);
yi = flip(cumsum(flip(y))) - y/2;
yi(1) = depth - f*h(1)/2;
vi = v./sqrt(depth)*sqrt(yi);
Fs = 374*h.*d.*cosd(phi).*vi.^2;
Fs(1) = Fs(1)*f;
end
chainpara.m
function [lc,mc,dc] = chainpara(id)
rhoFe = 7.9e3
rho = [3.2 7.0 12.5 19.5 28.12]
Lc = [78 105 120 150 180]*1e-3;
lc = Lc(id);mc = rho(id).*lc;
dc = 2*sqrt(mc/rhoFe/lc/pi);
end
solvequileq.m
function phi = solvequiled(Fb,Fw,Fs,M,f)
g = 9.81;N = length(Fb);%N杆的个数
rho = 1.025e3;rhoFe = 7.9e3;
[theta,phi,F] = deal(zeros(1,N));%角度初始化为0
for i = 1:N-1
fx = F(i)*sind(theta(i)) + Fs(i);%水平方向的力
if i == 1 ;fx = fx + Fw;end%浮标还受风的作用力
fy = Fb(i) + F(i)*cosd(theta(i));
if i==6;fy = fy - M*g + rho*(M/rhoFe)*g;end%第六个杆还受重物球的力
F(i+1) = hypot(fx,fy);%sqrt(fx.^2+fy.^2)合力
theta(i+1) = min(acosd(fy/F(i+1)),90);
end
phi = atan2d(F.*sind(theta)+Fs/2,F.*cosd(theta)+Fb/2);%力矩平衡
phi(phi>90) = 90;
phi(1) = atan2d(Fs(1)*f/2+Fw(1)*(f+(1-f)/2),Fb(1)*f/2);
plotmoor.m
function plotmoor(x,y,h,phi,Lc,cId,w,v,M,N,depth,xbuoy,f)
N = length(x);
%浮标
x1 = x(N) - cosd(phi(1));
x2 = x(N) + cosd(phi(1));
x3 = x(N-1) - cosd(phi(1));
x4 = x(N-1) - cosd(phi(1));
y1 = y(N) + sind(phi(1));
y2 = y(N) - sind(phi(1));
y3 = y(N-1) - sind(phi(1));
y4 = y(N-1) + sind(phi(1));
fill([x1 x2 x3 x4],[y1 y2 y3 y4]);
hole on
text(x(N-1) + 1.3,y(1:N-6,sprintf('2.1f m',f*h(1)));
%画锚链
plot(x(1:N-6),y(1:N-6),'m')
%画重物球
plot([x(N-6) x(N-6)],[y(N-6) y(N-6) - 4],'r')
plot(x(N-6),y(N-6) - 4,'.r','MarkerSize',50)
%标记风俗，水速
quiver(x(N) - 8,depth+1,2,0,'MaxHeadSize',0.5,'color','r');
text(x(N) - 5.5,depth+1,sprintf('%3.1fm/s',w),'color','r');
