求二元函数

z=0.2323*x^2-0.2866^2+2*(-0.5406)*a0^2+1.0203*a0^2*x^2/((x^2+y^2)^0.5*tanh(2*(x^2+y^2)^0.5)-x^2*(0.5733-u0)^2)

的最大值，变量x和y都是在0.2附近。求z的最大值，以及x,y的取值。

先用diff命令求z关于x,y的偏导数得到q(1)和q(2)两个方程，

再用fsolve求解方程组，想得到各驻点的坐标，但一直不成功，只得到x=0,y=0;

程序如下：

function q=myfun(p)

x=p(1);

y=p(2);

q(1)=(1046382468744793*x)/2251799813685248 + (2940675037869895*x)/(144115188075855872*(tanh(2*(x^2 + y^2)^(1/2))*(x^2 + y^2)^(1/2) - (2960289404442881*x^2)/9007199254740992)) + (2940675037869895*x^2*((2960289404442881*x)/4503599627370496 + 2*x*(tanh(2*(x^2 + y^2)^(1/2))^2 - 1) - (x*tanh(2*(x^2 + y^2)^(1/2)))/(x^2 + y^2)^(1/2)))/(288230376151711744*(tanh(2*(x^2 + y^2)^(1/2))*(x^2 + y^2)^(1/2) - (2960289404442881*x^2)/9007199254740992)^2);

q(2)=(2940675037869895*x^2*(2*y*(tanh(2*(x^2 + y^2)^(1/2))^2 - 1) - (y*tanh(2*(x^2 + y^2)^(1/2)))/(x^2 + y^2)^(1/2)))/(288230376151711744*(tanh(2*(x^2 + y^2)^(1/2))*(x^2 + y^2)^(1/2) - (2960289404442881*x^2)/9007199254740992)^2) - (1290927876507717*y)/2251799813685248；

main    % 主程序

x0 = [0.2;0.2];   % 初值

options = optimoptions('fsolve','Display','iter');

[x,fval] = fsolve(@myfun,x0,options),

谢谢。