1) You really should have simplified your conditions first by dividing out 2 and 4 before differentiating it. This does not change the result, but it is needless complication.
2) $x$, $y$ and $z$ are positive numbers, so we have to take into account the behaviour on the border when one or more of them are zero.
The situation is symmetric, so say $z=0$, but then you have $x+y=50$ and $xy= 750$ which gives (by the solution formula of the quadratic equation) that $x$ and $y$ are not real numbers.
This means that subject to the two conditions there is actually no border and every extrema can be found by using Lagrange multipliers.
3)
You write:
$x(y−z)/2=λ(y−z)$, so that $λ=x/2$
But what you get, for each of the three similar equations, is:
Either (a) $y=z$ or (b) $\lambda=x/2$.
Either (a) $z=x$ or (b) $\lambda=y/2$.
Either (a) $x=y$ or (b) $\lambda=z/2$.
Now, you would have to choose (a) or (b) for each of the three equations which gives 8 cases (aaa,aab,aba,abb,baa,bab,bba,bbb), but it is simpler than that.
If you choose at least once condition (a), then you have two equal variables. If you choose at least twice condition (b), then you have two variables equal to $2\lambda$, so you also have two equal variables.
So, whatever you do, two variables will be equal.
But the problem statement is symmetric in $x$, $y$ and $z$ (which means that if you interchange the three variables, the problem statement does not change, for example $yzx=xyz$, so the volume you are looking for will not depend on the order of the variables).
Therefore, since you know that two variables are equal and you know that it doesn't really matter for the end result which ones are equal, it suffices to investigate $x=y$:
(At this point, we still might have $x=y=z$, but we don't know that. What we do know, is that the original conditions are true, so we use them first, and it turns out that they already fix all the values.)
The original conditions give:
$2x+z=50$ and $x^2+2xz=750$.
Substituting $z=50-2x$ in the second equation gives a quadratic equation $3x^2-100x+750=0$.
This gives $x=\frac{50\pm 5\sqrt{10}}{3}$.
You find two solutions (up to symmetry):
$x=y=\frac53 (10 + \sqrt{10})$, $z=\frac{10}3(5-\sqrt{10})$ giving the volume
$\frac{2500}{27}(35-\sqrt{10})\approx 2948$
and
$x=y=\frac 53 (10- \sqrt{10})$, $z=\frac{10}3(5+\sqrt{10})$ giving the volume
$\frac{2500}{27}(35+\sqrt{10})\approx 3534$.
Since these are the only two candidates for extrema, the first one is the minimum value and the second one is the maximum value.
Suppose you want to maximize $z=f(x,y)$ subject to the constraint $g(x,y)=c$. You've used the method of Lagrange multipliers to have found the maximum $M$ and along the way have computed the Lagrange multiplier $\lambda$. Then $\lambda={dM\over dc}$, i.e. $\lambda$ is the rate of change of the maximum value with respect to $c$.
Said another way, you can think of $\lambda$ as approximately the change in $M$ that results from a one unit change in $c$.
Elaborating further: Optimizing $f(x,y)$ subject to $g(x,y)=c$ via Lagrange multipliers leads to $\nabla f=\lambda g$. Let $L(x,y;\lambda):=f(x,y)+\lambda(c-g(x,y))$. Then the constrained optimization problem can be cast as $\nabla L=0$.
From this perspective, ${\partial L\over \partial c}=\lambda$, i.e. $\lambda$ is the rate of change of the quantity being optimized, $L$, with respect to the constraint value, $c$.
Best Answer
Picking up from a little bit before where you left off, we can use the constraint to simplify the system of equations to $$y=\frac\lambda L x,x=\frac\lambda Ly$$ with solutions $\lambda=L$, $x=\pm y$. The variables $x$ and $y$ represent physical lengths, so the only solution we’re interested in is $x=y\gt0$, from which $x=y=L/\sqrt2$. This isn’t too surprising: the symmetry of both the objective function and constraint already told us that something interesting was likely to happen at $x=y$.
To determine the nature of these critical points, we examine second derivatives just as we might in elementary calculus. Here, we need to look at the behavior along the constraint curve, which we can do with the Hessian of the objective function. This is a quadratic form obtained by differentiating $f$ twice that we apply to the tangent to the constraint curve at the point of interest. $\nabla G=(2x,2y)$ and the tangent to a level curve is orthogonal to the gradient, so at $(x,y)=(L/\sqrt2,L/\sqrt2)$ a tangent vector is $(-L\sqrt2,L\sqrt2)$. We can use any vector that points in the same direction for this, so for simplicity we’ll use $(-1,1)$. Applying the Hessian to this vector we get $$H_f(-1,1)=\pmatrix{-1&1}\pmatrix{0&1\\1&0}\pmatrix{-1\\1}=-2\lt0$$ so there’s a local maximum along the curve at $x=y$.