2016年全国数模A题的问题1的MATLAB代码

代码如下：

%2016.9.11 系泊系统设计 
% r,H,v,V,beta,rou,L 
% r 是钢球质量比（以 1200kg 为基准），即钢球质量=1200kg*r，调整 r 即可调整钢球质量。H 为水深，v 为风速，  V为水速，beta 为风速水速之间的夹角，rou 是锚链密度，L 是锚链的长度 
%phi0,phi1,phi2,phi3,phi4,h,theta,x 
%phi0 是钢桶的倾斜角，phi1 至 phi4 为从下至上各节钢管的倾斜角，h 为吃水深度，theta为锚链末端与锚的连接 点处切线方向与海床夹角，x 为游动的最大范围。 
r=1; 
H=18; 
v=12; 
V=0; 
beta=0; 
rou=7; 
L=22.05; 
%考虑到锚链可能并不是全部悬空，有可能部分堆放在海床上，因此先求解不堆放情形，至于堆放情形，在下面 if 中考虑。 
%fzero函数是用于求解单变量非线性方程的，它需要提供一个初始值，并通过迭代逼近方程的根。fzero函数的使用方法比较简单，但是对于复杂的方程可能需要多次尝试不同的初始值才能得到正确的结果。

f=fzero(@(h)(h+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./...(-10248.12+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./...(-10073.98+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/...(-10073.98+78.28+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/...(-10073.98+78.28*2+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/...(-10073.98+78.28*3+31649.66.*h).^2)+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).*...(sqrt(1+(-22143.12+11760-11760*r*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2)...-sqrt(1+(-23655.75+7.*9.8.*22.05+11760-11760*r*0.87-rou.*L.*9.8.*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+...2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2))./(rou.*9.8.*0.87)-H),0.9); 
h=f; 

phi0=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10248.12+31649.66.*h).^2)); 
phi1=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./(-10073.98+31649.66.*h).^2)); 
phi2=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28+31649.66.*h).^2)); 
phi3=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28*2+31649.66.*h).^2)); 
phi4=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.Z^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28*3+31649.66.*h).^2)); 
theta=asec(sqrt(1+(-23655.75+7.*9.8.*22.05-rou.*L.*9.8.*0.87+11760-11760*r*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2)); 
alpha=asec(sqrt(1+(-22143.12+11760-11760*r*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2)); 
x=sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta))./(rou.*9.8.*0.87).*log((sec(alpha)+tan(alpha))/(sec(theta)+tan(theta)))+sin(phi0)+sin(phi1)+sin(phi2)+sin(phi3)+sin(phi4); 

%下面是判定锚链是否堆放在海床上，并进行此种情形的求解。 
outFywrong=(-22143.12+11760-11760*r*0.87+31649.66*h)%输出错误情形下 Fy 的值 
outGrou=rou.*9.8.*0.87.*L%输出错误情形下锚链的重力 
if sign(rou.*9.8.*0.87.*L-(-22143.12+11760-11760*r*0.87+31649.66*h))==1 

f=fzero(@(h)(h+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./...(-10248.12+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./...(-10073.98+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/...(-10073.98+78.28+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta))^2/...(-10073.98+78.28*2+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/...(-10073.98+78.28*3+31649.66.*h).^2)+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).*...(sqrt(1+(-22143.12+11760-11760*r*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2)...-1)./(rou.*9.8.*0.87)-H),0.7); 

h=f; 
phi0=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10248.12+31649.66.*h).^2)); 
phi1=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./(-10073.98+31649.66.*h).^2)); 
phi2=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28+31649.66.*h).^2)); 
phi3=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28*2+31649.66.*h).^2)); 
phi4=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28*3+31649.66.*h).^2)); 
theta=0; 
alpha=asec(sqrt(1+(-22143.12+11760-11760*r*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2)); 
x=sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta))./(rou.*9.8.*0.87).*log((sec(alpha)+tan(alpha)))+sin(phi0)+sin(phi1)+sin(phi2)+sin(phi3)+sin(phi4)+L-(31649.66*h-22143.12+11760-11760.*r.*0.87)./(rou.*9.8.*0.87); 
end 
outphi0=phi0.*180./pi%钢桶倾角 
outphi1=phi1.*180./pi%钢管倾角 
outphi2=phi2.*180./pi 
outphi3=phi3.*180./pi 
outphi4=phi4.*180./pi 
h 
outalpha=alpha.*180./pi 
outtheta=theta.*180./pi 
outx=x 
outk=sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta))./(rou.*9.8.*0.87)%锚链系数 k 
outFy=(-22143.12+11760-11760*r*0.87+31649.66*h) 
outFx=sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h))

放到MATLAB里运行，结果跟我说：

为啥呀？

有用望采纳：
我看了以下，应该是原来的代码中，fzero函数是用...加回车进行分割的，你复制的时候可能把回车给复制掉了。每个...后都加上回车就可以运行了。改成下面这样就可以了：

%2016.9.11 系泊系统设计 
% r,H,v,V,beta,rou,L 
% r 是钢球质量比（以 1200kg 为基准），即钢球质量=1200kg*r，调整 r 即可调整钢球质量。H 为水深，v 为风速，  V为水速，beta 为风速水速之间的夹角，rou 是锚链密度，L 是锚链的长度 
%phi0,phi1,phi2,phi3,phi4,h,theta,x 
%phi0 是钢桶的倾斜角，phi1 至 phi4 为从下至上各节钢管的倾斜角，h 为吃水深度，theta为锚链末端与锚的连接 点处切线方向与海床夹角，x 为游动的最大范围。 
r=1; 
H=18; 
v=12; 
V=0; 
beta=0; 
rou=7; 
L=22.05; 
%考虑到锚链可能并不是全部悬空，有可能部分堆放在海床上，因此先求解不堆放情形，至于堆放情形，在下面 if 中考虑。 
%fzero函数是用于求解单变量非线性方程的，它需要提供一个初始值，并通过迭代逼近方程的根。fzero函数的使用方法比较简单，但是对于复杂的方程可能需要多次尝试不同的初始值才能得到正确的结果。
 
f=fzero(@(h)(h+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./...
    (-10248.12+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./...
    (-10073.98+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/...
    (-10073.98+78.28+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/...
    (-10073.98+78.28*2+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/...
    (-10073.98+78.28*3+31649.66.*h).^2)+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).*...
    (sqrt(1+(-22143.12+11760-11760*r*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2)...
    -sqrt(1+(-23655.75+7.*9.8.*22.05+11760-11760*r*0.87-rou.*L.*9.8.*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+...
    2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2))./(rou.*9.8.*0.87)-H),0.9); 
h=f; 
 
phi0=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10248.12+31649.66.*h).^2)); 
phi1=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./(-10073.98+31649.66.*h).^2)); 
phi2=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28+31649.66.*h).^2)); 
phi3=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28*2+31649.66.*h).^2)); 
phi4=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28*3+31649.66.*h).^2)); 
theta=asec(sqrt(1+(-23655.75+7.*9.8.*22.05-rou.*L.*9.8.*0.87+11760-11760*r*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2)); 
alpha=asec(sqrt(1+(-22143.12+11760-11760*r*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2)); 
x=sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta))./(rou.*9.8.*0.87).*log((sec(alpha)+tan(alpha))/(sec(theta)+tan(theta)))+sin(phi0)+sin(phi1)+sin(phi2)+sin(phi3)+sin(phi4); 
 
%下面是判定锚链是否堆放在海床上，并进行此种情形的求解。 
outFywrong=(-22143.12+11760-11760*r*0.87+31649.66*h)%输出错误情形下 Fy 的值 
outGrou=rou.*9.8.*0.87.*L%输出错误情形下锚链的重力 
if sign(rou.*9.8.*0.87.*L-(-22143.12+11760-11760*r*0.87+31649.66*h))==1 
 
f=fzero(@(h)(h+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./...
    (-10248.12+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./...
    (-10073.98+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/...
    (-10073.98+78.28+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta))^2/...
    (-10073.98+78.28*2+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/...
    (-10073.98+78.28*3+31649.66.*h).^2)+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).*...
    (sqrt(1+(-22143.12+11760-11760*r*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2)...
    -1)./(rou.*9.8.*0.87)-H),0.7); 
 
h=f; 
phi0=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10248.12+31649.66.*h).^2)); 
phi1=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./(-10073.98+31649.66.*h).^2)); 
phi2=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28+31649.66.*h).^2)); 
phi3=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28*2+31649.66.*h).^2)); 
phi4=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28*3+31649.66.*h).^2)); 
theta=0; 
alpha=asec(sqrt(1+(-22143.12+11760-11760*r*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2)); 
x=sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta))./(rou.*9.8.*0.87).*log((sec(alpha)+tan(alpha)))+sin(phi0)+sin(phi1)+sin(phi2)+sin(phi3)+sin(phi4)+L-(31649.66*h-22143.12+11760-11760.*r.*0.87)./(rou.*9.8.*0.87); 
end 
outphi0=phi0.*180./pi%钢桶倾角 
outphi1=phi1.*180./pi%钢管倾角 
outphi2=phi2.*180./pi 
outphi3=phi3.*180./pi 
outphi4=phi4.*180./pi 
h 
outalpha=alpha.*180./pi 
outtheta=theta.*180./pi 
outx=x 
outk=sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta))./(rou.*9.8.*0.87)%锚链系数 k 
outFy=(-22143.12+11760-11760*r*0.87+31649.66*h) 
outFx=sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h))

MATLAB中使用=运算符是用于赋值操作的，而不是用于比较相等性的判断。因此，将=替换为==即可修复这个错误。

请将代码中的h=f修改为h==f，然后再次运行代码，看看是否能够解决该错误。

该代码的作用是解决2016年全国数模A题问题1中的系泊系统设计问题。

具体实现原理如下： 1. 首先定义了一些需要用到的常量和变量。 2. 使用fzero函数求解非线性方程，其中方程是关于h的，通过迭代逼近得到方程的解。这个方程是用来计算不堆放情况下的吃水深度。 3. 接下来进行绘制图表的操作，使用cumsum函数求解累计和，并绘制出x和y的图像。其中，x是吃水深度的累计和，y是水深的累计和。 4. 最后输出一些计算结果。

代码具体实现了解决2016年全国数模A题问题1中的系泊系统设计问题，通过计算得到了不堆放情况下的吃水深度以及绘制了相应的图表。