Note first that since $L(x,y,z,k) = −xyz−k(xy+2z(x+y)−50) = L(y,x,z,k)$, the problem is symmetric in $x,y$.
Let's proceed with brutal arithmetic. We have 4 equations in 4 unknowns:
$$
\begin{split}
0 &= -L_x = yz + k(y + 2z)\\
0 &= -L_y = xz + k(x + 2z)\\
0 &= -L_z = xy + 2k(x + y)\\
50 &= xy+2z(x+y) \quad (\text{from } c=0)
\end{split}
$$
Claim 1. $x \neq 0$ and $y \neq 0$ and $x + y \neq 0$.
Proof.
Note first that if $x=0$, then the second constraint implies $k=0$ or $z=0$. If $x=k=0$, the first constraint implies $yz=0$ and the last constraint implies $2yz=50$, which is impossible.
Similarly, if $x=z=0$, the first constraint implies $ky = 0$. As we proved, $k=0$ is impossible, so must be $x=y=z=0$, which contradicts the last constraint.
To establish the last claim, note that if $x+y =0$, the third constraint yields $xy=0$ and the last one yields $xy=50$, which is a contradiction.
QED
Now, solving the last constraint for $z$ and the $L_z$ constraint for $k$, we get
$$k = \frac{-xy}{2(x+y)} \text{ and } z = \frac{50-xy}{2(x+y)}.$$
Both of these are well-defined since $x\neq0$ and $y \neq 0$. Plug both of these into the first constraint, getting
$$0 = \frac{y(50-xy)}{2(x+y)}
- \frac{xy^2}{2(x+y)}
- \frac{2xy(50-xy)}{2^2(x+y)^2}$$
and now multiply both sides by $2(x+y)^2 \neq 0$ to get
$$0 = y(50-xy)(x+y) - xy^2(x+y) - xy(50-xy).$$
Divide by $y \neq 0$ and bring the last term into the first term and divide by $y$ again:
$$
\begin{split}
0 &= (50-xy)(x+y) - xy(x+y) - x(50-xy)\\
0 &= y(50-xy) - xy(x+y) \\
0 &= 50-xy - x(x+y) = 50 - 2xy - x^2\\
50 &= x^2 + 2xy
\end{split}
$$
and because the problem is symmetric, the symmetric constraint must hold (if you want, you can derive it the same way from the second constraint):
$$
\begin{split}
x^2 + 2xy &= 50 \\
y^2 + 2xy &= 50
\end{split}
$$
Hence subtracting them yields $x^2 = y^2$, so $x = \pm y$, but $x \neq -y$ by Claim 1, so $x = y$. Now $50 = x^2 + 2xy = 3x^2$ implies $x = \pm \sqrt{50/3} = y$, and there are two solutions: $(\sqrt{50},\sqrt{50})$ and $(-\sqrt{50},-\sqrt{50})$.
Since $x,y$ are lengths, we know only one will make sense: $x = y = \sqrt{50/3}$. Now plug this back in to find $z$ and $k$ as desired.
Consider a point $p$ in the common domain $\Omega\subset{\mathbb R}^n$ of $f$ and the constraints $$g_k(x)=0\qquad(1\leq k\leq r)\ .\tag{1}$$ The gradients $\nabla g_k(p)$ define a subspace $U$ of allowed directions when walking away from $p$. In fact a direction $X$ is allowed only if it belongs to the tangent planes of all level surfaces $(1)$. This means that $X$ is perpendicular to all $\nabla g_k(p)$, or is a solution of the homogeneous system of equations
$$\nabla g_k(p)\cdot X=0\qquad(1\leq k\leq r)\ .\tag{2}$$
Now comes an important technical condition for the application of Lagrange's method: We have to assume that the $r$ gradients $\nabla g_k(p)$ are linearly independent, i.e. that $p$ is a regular point of the manifold defined by $(1)$. In this case the $\nabla g_k(p)$ span an $r$-dimensional subspace $V$, and the system $(2)$ has full rank. It follows that ${\rm dim}(U)= n-r$. Therefore we not only have $U\subset V^\perp$, but in fact
$$U=V^\perp\ .$$
When $\nabla f(p)\cdot X\ne0$ for some allowed direction $X$ then the function $f$ is not conditionally stationary at $p$. For a constrained local extremum of $f$ at $p$ we therefore need $$\nabla f(p)\cdot X=0$$
for all directions $X\in U$, in other words: It is necessary that $$\nabla f(p)\in U^\perp=V\ .\tag{3}$$
When $\nabla f(p)\in V=\langle\nabla g_1(p),\ldots,\nabla g_r(p)\rangle$ then there are numbers $\lambda_k$ $\>(1\leq k\leq r)$ such that
$$\nabla f(p)=\sum_{k=1}^r \lambda_k\>\nabla g_k(p)\ .\tag{4}$$
Solving $(4)$ (with $x$ in place of $p$) together with $(1)$ will bring all regular constrained extrema of $f$ to the fore.
Best Answer
A problem could be
$\texttt{max} \ f(x,y)=-(x-1)^2-(y-1)^2$ under the constraints $x+y\leq 1$ and $x,y\geq 0.$
The lagrange function then is:
$L(x,y,\lambda )=-(x-1)^2-(y-1)^2+\lambda (1-x-y)$
The expression in the brackets of $\lambda ()$ has to be greater or equal to zero.
The KKT conditions are:
$\frac{\partial L}{\partial x}=-2(x-1)-\lambda\leq 0 \quad (1), \quad \frac{\partial L}{\partial y}=-2(y-1)-\lambda \leq 0 \quad (2)$
$\frac{\partial L}{\partial \lambda}=1-x-y\leq 0\quad (3), \quad x\cdot \frac{\partial L}{\partial x}=-x\left( 2(x-1)+\lambda\right)= 0\quad (4)$
$y\cdot \frac{\partial L}{\partial y}=-y(2(y-1)+\lambda) = 0 \quad (5), \quad \lambda\cdot \frac{\partial L}{\partial \lambda}=\lambda(1-x-y)=0 \quad (6), \quad x,y,\lambda \geq 0 \quad (7)$
Now you check the two cases $\lambda=0$ and $\lambda \neq 0 $
For (4) and (5) you would have 4 possible solutions $(x,y,\lambda)$:$(0,0,0),(1,0,0),(0,1,0),(1,1,0) $
None of these solutions satisfies the conditions (1),(2) and (3) simultaneously.
Because of (6) we have $1-x-y=0$
If $x=0$, then $y=1$. Inserting the values in (5): $-1(0-\lambda)=\lambda=0$ Because of $\lambda \neq 0$ (case 2) we have a contradiction.
If $x=1$, then $y=0$. Inserting the values in (4): $-1(0-\lambda)=\lambda=0$ Because of $\lambda \neq 0$ (case 2) we have a contradiction.
We can conclude, that $x,y \neq 0$. Because of (4) and (5) we have the two equations.
$2x-2-\lambda=0$
$2y-2-\lambda=0$
Substracting the second equation by the first equation.
$2x-2y=0 \Rightarrow 2x=2y \Rightarrow x=y$
With (6) and $\lambda\neq 0$ we get $1-x-x=0 \Rightarrow 1=2x \Rightarrow x=\frac{1}{2}$ and $y=\frac{1}{2}$ By using (5) we can calculate $\lambda$.
$2\cdot \left(\frac{1}{2}-1\right)+\lambda=0\Rightarrow \lambda=1$