How can I solve this complex systems of equations?

 Declaring variables a1, a2, b1, and b2 (numeric representation of letters)
 a1 = 16;
 a2 = 1;
 b1 = 23;
 b2 = 8;
 p = 4*(a1 + a2 + b1 + b2);

 %Formula for landscape 
 %creating system variables
 syms x y z;

 f = (p/((3*x+a1).^2 + (3*y+b1).^2+(a1+b1)))...
    -(2*p/((3*x-a2).^2 + (3*y+b2).^2+(a2+b2)))...
    -(3*p/((3*x+a1).^2 + (3*y-b2).^2+(a1+b2)))...
    +(4*p/((3*x-a2).^2 + (3*y-b1).^2+(a2+b1)));

 %plotting the function from -10 to 10 in x&y domain
 ezsurf(f, [-10, 10]);

 %first partial derivatives w/respect to x & y

 fx = diff(f, x)==0;
 fy = diff(f, y)==0;
 [xcr, ycr] = solve (fx, fy); [xcr, ycr];

the script gets held up on the "solve" command and freezes. I have no idea where to go from here as I need to find the critical points in order to find the local min, max, and saddle points of the function. any help is much appreciated

Answer 1

[Star Strider]

This is likely the best you can hope for, considering that your derivative equations end up to be high-degree polynomials in both ‘x’ and ‘y’:

fx = simplify(diff(f, x),'Steps',20);
fy = simplify(diff(f, y), 'Steps',20);

fxf = matlabFunction(fx);
fyf = matlabFunction(fy);
fxyf = @(b) [fxf(b(1),b(2)); fyf(b(1),b(2))];
B = fsolve(fxyf, [0 0]);

You will need to loop through appropriate initial values for ‘b(1)=x’ and ‘b(2)=y’. Use the uniquetol function to eliminate duplicate values for ‘b’ with 'ByRows',true.

[Patrick Wheeler]

Thanks! So, I understand the process of simplifying the expression up until the declaration of "fxyf". What exactly is happening here?

Also, could you explain what you mean by looping through b(1) and b(2) using uniquetol? The resultant of the fsolve only gives two different values, so there aren't any duplicates to eliminate.

Additionally, based graph, I would've expected at least four values(two local min, two local max, and saddle points). The two that are given seem to indicate the saddle points...what am I missing?

Thanks again for taking the time to address this question, it means a lot!

[Star Strider]

My pleasure!

My apologies for not documenting my code. The ‘fxyf’ function combines the two anonymous function calls to ‘fx’ and ‘fy’ into one function to give to fsolve, that finds the solutions closest to the initial parameter estimates (here, the ‘parameters’ are ‘x’ and ‘y’).

You have to loop through the initial parameter estimates around the known maxima and minima to get fsolve to identify them. I would use a for loop, with perhaps 5 steps between ±5 for the initial parameter estimates (so [-5, -5] to [5 5], most easily using two nested loops, one loop for each variable), saving each set of parameters estimated by fsolve as a separate row in your output array, then use uniquetol to eliminate the duplicates.

Example:

k = 0;
for b1 = -5 : 1 : 5
    for b2 = -5 : 1 : 5
        k = k + 1;
        B(k,:) = fsolve(fxyf, [b1 b2]);
    end
end
crit_pts = uniquetol(B, 1E-1, 'ByRows',true);

You will likely have to experiment with it to get the result you want.

Answer 2

[Walter Roberson]

Solution attached (it was too long to post.)

Afterwards, xy will be the symbolic solutions and xyd will be the numeric equivalents.

xyd will have x values in column 1, y values in column 2. 6 of the rows are purely real-valued and the rest are imaginary valued.

The symbolic solutions are pretty much useless, being the root of a 36 degree polynomial.