我有一个13方程13未知数的方程组需要求解，使用的是fsolve函数，设定初值之后每次计算的结果都不同且相差比较大，精度不高。

希望各位大神能帮我修改一下计算程序以提高计算精度。

下面是自定义函数的代码：

function q=myfun(p)

t3=13;t5=58;t8=94;t9=80;h1=2666.1;h2=305.6;h3=2405.9;D=3.643;%已知系数

%下面是需要求解的方程组，一共13个方程，13个未知数。

x4=p(1); x5=p(2); x6=p(3); x8=p(4); x9=p(5); t4=p(6); t6=p(7); h4=p(8); h5=p(9); h6=p(10); h8=p(11); h9=p(12); G=p(13);

q(1)=(G-D)*(h8-h9)-G*(h6-h5);

q(2)=(G-D)*x8-G*x6;

q(3)=D*h1+(G-D)*h8-G*h6+D*(h3-h2)-20000;

q(4)=(h8-h9)/(h8-h4)-0.65;

q(5)=G*t4-D*t3-(G-D)*t9;

q(6)=G*x5-(G-D)*x9;

q(7)=x8-x9;

q(8)=x5-x6;

q(9)=-2024.33+163.309*x4-4.88161*x4^2+(6.302948e-2)*x4^3-(2.913705e-4)*x4^4+t4*(18.2829-1.1692*x4+(3.24804e-2)*x4^2-(4.0342e-4)*x4^3+(1.85206e-6)*x4^4)+t4^2*(-(3.7008e-2)+(2.88777e-3)*x4-(8.1313e-5)*x4^2+(9.91166e-7)*x4^3-(4.44412e-9)*x4^4)-h4;

q(10)=-2024.33+163.309*x5-4.88161*x5^2+(6.302948e-2)*x5^3-(2.913705e-4)*x5^4+t5*(18.2829-1.1692*x5+(3.24804e-2)*x5^2-(4.0342e-4)*x5^3+(1.85206e-6)*x5^4)+t5^2*(-(3.7008e-2)+(2.88777e-3)*x5-(8.1313e-5)*x5^2+(9.91166e-7)*x5^3-(4.44412e-9)*x5^4)-h5;

q(11)=-2024.33+163.309*x6-4.88161*x6^2+(6.302948e-2)*x6^3-(2.913705e-4)*x6^4+t6*(18.2829-1.1692*x6+(3.24804e-2)*x6^2-(4.0342e-4)*x6^3+(1.85206e-6)*x6^4)+t6^2*(-(3.7008e-2)+(2.88777e-3)*x6-(8.1313e-5)*x6^2+(9.91166e-7)*x6^3-(4.44412e-9)*x6^4)-h6;

q(12)=-2024.33+163.309*x8-4.88161*x8^2+(6.302948e-2)*x8^3-(2.913705e-4)*x8^4+t8*(18.2829-1.1692*x8+(3.24804e-2)*x8^2-(4.0342e-4)*x8^3+(1.85206e-6)*x8^4)+t8^2*(-(3.7008e-2)+(2.88777e-3)*x8-(8.1313e-5)*x8^2+(9.91166e-7)*x8^3-(4.44412e-9)*x8^4)-h8;

q(13)=-2024.33+163.309*x9-4.88161*x9^2+(6.302948e-2)*x9^3-(2.913705e-4)*x9^4+t9*(18.2829-1.1692*x9+(3.24804e-2)*x9^2-(4.0342e-4)*x9^3+(1.85206e-6)*x9^4)+t9^2*(-(3.7008e-2)+(2.88777e-3)*x9-(8.1313e-5)*x9^2+(9.91166e-7)*x9^3-(4.44412e-9)*x4^4)-h9;

下面是主程序代码：

clear all

clc

[x,fval,]=fsolve(@myfun,rand(1,13)*100,optimset('MaxIter',1e5,'MaxFunEvals',1e6))