I want to find the maximum of the following function, but get a wrong result. What am I doing wrong? Clear[x, y, xx, yy] peaks[x_, y_] := 3*(1 - x)^2.*Exp[-(x^2) - (y + 1)^2] - 10*(x/5 - x^3 - y^5) Exp[-x^2 - y^2] - 1/3 Exp[-(x + 1)^2 - y^2]; N[Maximize[peaks[xx, yy], {xx, yy}]] NMaximize[{peaks[x, y], x >= -2, x = -2, y peaks[0, 2] peaks[1.28568, -0.00484756] When you look at these graphs, you realize that the results are strange: peaks[x_, y_] := 3*(1 - x)^2*Exp[-(x^2) - (y + 1)^2] - 10*(x/5 - x^3 - y^5) Exp[-x^2 - y^2] - 1/3 Exp[-(x + 1)^2 - y^2]; Plot3D[peaks[x, y], {x, -2, 2}, {y, -2, 3}, ColorFunction -> "DarkRainbow", AxesLabel -> Automatic] ContourPlot[peaks[x, y], {x, -2, 2}, {y, -2, 3}, ColorFunction -> "DarkRainbow", ContourLines -> False, Contours -> 25] Answer You should use FindMaximum : FindMaximum[peaks[xx, yy], {xx, yy}] {8.106, {xx -> -0.009, yy -> 1.581}} plot = Plot3D[peaks[x, ...