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,