Maxima: unable to solve Lagrangian I think this is may be more of a technical issue than a mathematical issue. I don't know...
I have the following function I want to maximize:
$$\arg\max_{x,y}\ U(x,y) = \frac{1-e^{-x}}{2} + \frac{1-e^{-y}}{2}$$
subject to $x+y \leq 50$ and $x,y\geq0$.
(My actual poblem is more complicated than this of course, but take this as a MWE please. :))
We can easily see this has a maximum. Using Maxima:
U(x,y) := (1-exp(-x))/2 + (1-exp(-y))/2;
wxplot2d([U(x,50-x)], [x,0,50]);


To solve this using Lagrangian optimzation along the $x+y=50$ frontier, I would do:
$$\begin{cases}
\dfrac{\partial U}{\partial x} = \lambda \\[4pt]
\dfrac{\partial U}{\partial y} = \lambda \\[4pt]
x+y = 50
\end{cases}
$$
We can easily solve this. $\frac{\partial U}{\partial x} = \frac{e^{-x}}{2}=\lambda \iff \frac{e^{-x}}{2} = \frac{e^{-50+x}}{2} \iff x=25$, as could be seen in the graph anyway.
But now I try to do this in Maxima:
U(x,y) := (1-exp(-x))/2 + (1-exp(-y))/2;
solve([diff(U(x,y),x)=lambda, diff(U(x,y),y)=lambda, x+y=50], [x,y,lambda]);

And it says there is no solution:
[]

What gives? I haven't done anything in optimization or Maxima for a long time... But I cannot for the life of me discover what the problem is...
EDIT: I would prefer to know where my Maxima usage is wrong (or my mathematics), but I will approve any answers that suggest me to use another symbolic solver.
 A: Your solution is correct. It is also the same in Mathematica. One can still use "Reduce" in Mathematica as follows:
Reduce[{E^-x/2 - a == 0, E^-y/2 - a == 0, x + y == 50}, {x, y, a}]

with the following answer:
C[1] \[Element] Integers && (x == 25 + I \[Pi] + 2 I \[Pi] C[1] || x == 25 + 2 I \[Pi] C[1]) && y == 50 - x && a == E^-x/2

Since the problem is simple one can solve this problem without using any Lagrangian multiplier as well.
The function $U$ is increasing both in $x$, $y$ and in both. So this means for the maximization $x$ and $y$ should be as large as possible. So this means the constraint $x+y\leq 50$ should hold with equality, that is $x+y=50$. Then, one can simply insert $x:=-y+50$ in $U$, take the derivative w.r.t. $y$ and see that $y$ must be $25$. Accordingly $x=25$.
You may think about asking questions about programming in mathematica.stackexchange
A: In Maple:
U:= (1-exp(-x))/2 + (1-exp(-y))/2;
solve({diff(U,x)=lambda, diff(U,y)=lambda, x+y=50},{lambda,x,y});

$$ \left\{ \lambda=\frac{1}{2 {\rm e}^{25}},x=25,y=25
 \right\} 
 $$
A: For future reference, I was able to solve this using Sympy, though I had to get rid of the $\lambda$... (If someone knows why please let me know.)
from sympy import *
x, y = symbols('x y')

U = (1-exp(-x))/2 + (1-exp(-y))/2
dU_dx = diff(U, x)
dU_dy = diff(U, y)

print solve((Eq(dU_dx,dU_dy), Eq(x+y,50)), (x, y))

