Analytic solution for the maxima of a bivariate

Question

I found the maximum of the function $f(x) = \frac{x e^{-x}}{1+1/k-e^{-x}}$ by reducing the first order necessary condition to $ke^{-x}+(1+k)(x-1)=0$, and from there the solution obtained with a special (Lambert-W) function at $\;x = \mathbb{W}_0\left(-\frac{k}{e(k+1)}\right)+1$ (thanks to Solve[] function in Mathematica).
Now I have a more complicated, bi-variate version of the function to deal with: $$ f(x,y) = \frac{(k+1)e^{-x}(xe^{x+y}+y)}{(k+2)(e^{x+y}-1)+e^y} $$ where $k \in \mathbf{N}$ is a constant, and both $x,\, y \in \mathbf{R_+}$. Maximum of the function appears in the plot to remain for small values of $(x, y)$. So far I haven't had any luck in solving $\frac{\partial f}{\partial x} = 0,\; \frac{\partial f}{\partial y} = 0$ as both involve exponential forms (including Mathematica's Solve[]).
I need an analytic solution, even if it uses special function(s) and/or is a close approximation. Any ideas would be highly appreciated.