Hi,
my machine also doesnt stop being busy, when trying to solve this problem symbollically, so here is an numeric approach:
syms f(x,y)
fun = f == (100/((3*x+13)^2+(3*y+2)^2+(13+2))) - ((2*100)/((3*x+5)^2+(3*y-5)^2+10)-((3*25)/((3*x+13)^2+(3*y-5)^2+18))+((4*100)/(3*x-5)^2+(3*y-2)^2+7));
fsurf(rhs(fun),[-10,10])
colormap(jet)
hold on
% define the function f as function handle (same with the derivates at the
% function below. x0 is a guess from looking at the function plot.
func = @(x)-(1.0./(x(1).*3.0-5.0).^2.*-4.0e2-(x(2).*3.0-2.0).^2-2.0e2./((x(1).*3.0+5.0).^2+(x(2).*3.0-5.0).^2+1.0e1)+1.0e2./((x(1).*3.0+1.3e1).^2+(x(2).*3.0+2.0).^2+1.5e1)+7.5e1./((x(1).*3.0+1.3e1).^2+(x(2).*3.0-5.0).^2+1.8e1)-7.0);
sol = fmincon(func,[-5 0.5],[],[],[],[],[-10 -10], [10 10], @derivates)
% show the solution on the plot
max_value = -func(sol);
scatter3(sol(1), sol(2), max_value, 'og', 'LineWidth',3)
hold off
% Make sure that the derivates of f with respect to x and y
% are equal to zero
function [c,ceq] = derivates(x)
c = [];
ceq = [(x(1).*1.8e1+3.0e1).*1.0./((x(1).*3.0+5.0).^2+(x(2).*3.0-5.0).^2+1.0e1).^2.*2.0e2-(x(1).*1.8e1+7.8e1).*1.0./((x(1).*3.0+1.3e1).^2+(x(2).*3.0+2.0).^2+1.5e1).^2.*1.0e2-(x(1).*1.8e1+7.8e1).*1.0./((x(1).*3.0+1.3e1).^2+(x(2).*3.0-5.0).^2+1.8e1).^2.*7.5e1+1.0./(x(1).*3.0-5.0).^3.*2.4e3, x(2).*-1.8e1-(x(2).*1.8e1+1.2e1).*1.0./((x(1).*3.0+1.3e1).^2+(x(2).*3.0+2.0).^2+1.5e1).^2.*1.0e2+(x(2).*1.8e1-3.0e1).*1.0./((x(1).*3.0+5.0).^2+(x(2).*3.0-5.0).^2+1.0e1).^2.*2.0e2-(x(2).*1.8e1-3.0e1).*1.0./((x(1).*3.0+1.3e1).^2+(x(2).*3.0-5.0).^2+1.8e1).^2.*7.5e1+1.2e1];
end
By using this i got:
sol =
-4.727519221259621 0.605820153518793
which gives this as result when calculating the function value:
>> max_value
max_value =
-4.266365719653375
or checked against the result fro your symbolic function:
>> check = double(subs(rhs(fun),[x y],sol))
check =
-4.266365719653376
As far as in the plot can be seen, this should be (one) correct solution to your function. Note that it could be a local extremum.
Best regards
Stephan
