How to find the symbolic maximum of a 3D function?

edit

You have two types of conditions:

• $a>0$, $b>0$ : those are assumptions on the parameters, not part of a condition on candidate solutions

• $x>0$, $y>0$, $y<=x$ : those are restrictions imposed to candidate solutions.

Therefore :

# Problem-specific setup
var("x,y,a,b", domain="real")
f(x,y)=a*(x*y-y^2)*e^(-b*x^2)*e^(b*x*y)
assume(a>0,b>0)
# Extra conditions
Conds=[x>0, y>0, y<=x]
# Derivatives : must be all null to have an extremum/saddle point
Ders=[f(*f.arguments()).diff(u) for u in f.arguments()]
# Second derivatives : must be all negative to have a maximum
SDers=[f(*f.arguments()).diff(u,2) for u in f.arguments()]
# Extrema/saddle points
Exts=solve(Ders,[x,y], to_poly_solve="force", solution_dict=True)
# Candidate solutions = maxima among extrema
CSols=[u for u in Exts if all([bool(s.subs(u)<=0) for s in SDers])]
# Solutions : maxima fulfilling extra conditions
Sols=[u for u in CSols if all([v.subs(u) for v in Conds])]

sage: Sols
[]

Note that

sage: CSols
[{x: 0, y: 0}]

The (only) maximum $f$ over $\mathbb{R}^2$ is excluded by your restrictions.

HTH,

EDIT : I've found an easier way yo show this. Consider your problem in polar coordinates, or, more lazily, by posing y=k*x (i. e. k=arctan(theta)).

sage: var("a, b, x, y, k", domain="real")
(a, b, x, y, k)
sage: f(x, y) = a * (x * y - y^2) * exp( -b * x^2) * exp(b * x * y )

We can get rid of the special cases of the axes by noting :

sage: f(x,0)
0
sage: f(0,y)
-a*y^2

The former is constant, the latter is maximal (and zero) in y=0.

Now, what is the maximum of f along the line of equation y=k*x ?

sage: Sx=f(x,k*x).diff(x).solve(x);Sx
[x == -sqrt(-1/(b*k - b)), x == sqrt(-1/(b*k - b)), x == 0]

(The third case has already been disposed of...). Now, what are the values of k maximizing f . Well, the derivatives of these maxima with respect to k :

sage: Dk=[f(x,k*x).subs(u).diff(k).simplify_full() for u in Sx[:2]];Dk
[a*e^(-k/(k - 1) + 1/(k - 1))/b, a*e^(-k/(k - 1) + 1/(k - 1))/b]

are constant :

sage: [u.log().expand_log().factor().exp().simplify_log() for u in Dk]
[a*e^(-1)/b, a*e^(-1)/b]

Therefore no bloody maximum, notwithstanding the numerical phantasms of Sage (or Mathematica, for that matter...).

HTH,