求大佬指教！！用fsolve函数求解非线性方程组

huangdongqiang

function equation()
global t1  L1 L2 L3  a1 a2 a3 E  I  LX LY d H t2 w1 w2
L1=15e-3;
L2=15e-3;
L3=2e-3;
LX=31.514e-3;
a1=10*pi/180;
a2=10*pi/180;
a3=10*pi/180;
t1=1.2e-3;
t2=1.2e-3;
LY=5.557e-3;
E=1.420e9;
H=6e-3;
w1=0;
w2=0;
I=(t1)^3*H/12;
N=1;
for d=0:0.0001:8e-3   
    %x0=[0 0 0 0 0 0 0 0 0 0 0 0 0]*2;
    x0=rand([14,1])*1;
    % 定义迭代初值
    %opt=optimset('Display','off','MaxFunEvals',1e20,'MaxIter',1e20,'TolFun',1e-20,'TolX',1e-20);
    %options = optimoptions('fsolve','Display','none','PlotFcn',@optimplotfirstorderopt);
    %x0=[0,1];
    options = optimoptions('fsolve','Display','iter');
    [x,fval]=fsolve(@myfun,x0,options);    % 采用fsolve函数解对应lambda下的方程组，结果保存在x里
    value(N)=6*x(13);   
    N=N+1;
end
d=0:0.0001:8e-3;
plot(d,value)   % 绘图


end


function F=myfun(x)
global t1  L1 L2 L3  a1 a2 a3 E  I  LX LY d H t2 w1 w2
% x (1)~x(13)分别对应13个未知数
f1=x(1);
p1=x(2);
m1=x(3);
thea1=x(4);
dx1=x(5);
dy1=x(6);
f2=x(7);
p2=x(8);
m2=x(9);
thea2=x(10);
dx2=x(11);
dy2=x(12);
F=x(13);
% 定义13个方程
eq1=12*dy1-6*thea1+p1*(1.2*dy1-0.1*thea1)-f1+p1^2*[-1/700 1/1400]*[dy1;thea1];
eq2=-6*dy1+4*thea1+p1*(-0.1*dy1+2/15*thea1)-m1+p1^2*[1/1400 -11/6300]*[dy1;thea1];
eq3=t1^2*p1/(12*L1^2)+[dy1 thea1]*[-0.6 0.05;0.05 -1/15]*[dy1;thea1]-dx1+p1*[dy1 thea1]*[1/700 -1/1400;-1/1400 11/6300]*[dy1;thea1];
eq4=[12 -6]*[dy2;thea2]+p2*[1.2 -0.1]*[dy2;thea2]-f2+p2^2*[-1/700 1/1400]*[dy2;thea2];
eq5=[-6 4]*[dy2;thea2]+p2*[-0.1 2/15]*[dy2;thea2]-m2+p2^2*[1/1400 -11/6300]*[dy2;thea2];
eq6=t2^2*p2/(12*L2^2)+[dy2 thea2]*[-0.6 0.05;0.05 -1/15]*[dy2;thea2]-dx2+p2*[dy2 thea2]*[1/700 -1/1400;-1/1400 11/6300]*[dy2;thea2];
eq7=((p1*E*I/L1^2)*cos(a1)-(f1*E*I/L1^2)*sin(a1))+((p2*E*I/L2^2)*cos(pi+a2)-(f2*E*I/L2^2)*sin(pi+a2));
eq8=((p1*E*I/L1^2)*sin(a1)-(f1*E*I/L1^2)*cos(a1))+((p2*E*I/L2^2)*cos(pi+a2)-(f2*E*I/L2^2)*sin(pi+a2));
eq9=(m1*E*I/L1)+(m2*E*I/L2)-(f1*E*I/L1^2)*L3*cos(a3-a1-thea1)+(p1*E*I/L1^2)*L3*sin(a3-a1-thea1);
eq10=LX-((L1+(dx1*L1))*cos(a1)-dy1*L1*sin(a1)+L3*cos(a3)-((L2+dx2*L2)*cos(a2+pi)-dy2*L2*sin(a2+pi)));
eq11=(LY-d)-((L1+dx1*L1)*sin(a1)+dy1*L1*cos(a1)+L3*sin(a3+thea1)-((L2+dx2*L2)*sin(a2+pi)+dy2*L2*cos(a2+pi)));
eq12=F-((p2*E*I/L2^2)*sin(a2+pi)+(f2*E*I/L2^2)*cos(a2+pi));
eq13=thea1-thea2;
% 返回方程组
F=[eq1;eq2;eq3;eq4;eq5;eq6;eq7;eq8;eq9;eq10;eq11;eq12;eq13];
end
这个结果是



这个波动的太大了 请问怎样调整可以减小波动？使曲线更加平滑，