Let $L$ be defined as follows
$$L = \{ (x,y,z) : x - y + 2z =4 , \; \; \; 2x+y-z = 1 \} $$
I want to find the point in $L$ closest to $(0,0,0)$
LEt $(x,y,z)$ be arbitrary point in $L$, the distance of this point from $0$ is
$$ F = D^2 = x^2 + y^2 + z^2 $$
So, I want to minimize this function subject to the constrains
$g_1 = x - y + 2z - 1 = 0$ and $g_2 = 2x + y - z - 1 = 0$.
Is this a correct approach?
You asked for an answer from a credible source. I used to be a professor of physics and mathematics, although that was a while ago now so I leave the judgment of credibility up to you.
You also asked for a solution with Lagrange multipliers. Here is one.
Construct the objective function $$F(x,y,z) = x^2+y^2+z^2 + \lambda(x-y+2z) + \mu(2x + y - z)\,.$$ Here $\lambda$ and $\mu$ are the Lagrange multipliers. We need two because there are two constraints to satisfy.
Take partial derivatives of $F$, set them equal to zero, and find the following results for $x$, $y$, and $z$:
$$x = -{\lambda \over 2} - \mu$$ $$y = {\lambda - \mu\over 2}$$ $$z = {\mu\over 2} - \lambda\,.$$
By the second derivative test, these choices will give a local minimum (regardless of the values of $\lambda$ and $\mu$).
Impose the two constraint equations using these expressions for $x$, $y$, and $z$. This leads to
$$6\lambda - \mu = -8$$ $$\lambda - 6\mu = 2\,.$$
Solve these equations to find $\lambda = -{10\over 7}$ and $\mu = -{4\over 7}$. These go into the expressions for $x$, $y$, and $z$ above to yield
$$x = {9\over 7}$$ $$y = -{3\over 7}$$ $$z = {8\over 7}\,.$$
And these lead to $$F = {22\over 7}\,.$$
You can avoid constrained optimization. Write the equation of a line in parametric form. (Solve the linear non-homog. system and get the infinitely many solutions expressed in terms of a single free variable). So a point on the line will be of the form $(a_1+b_1t, a_2+b_2t, a_3+b_3t)$ Now the square of the distance of this point to origin (or any other point) is a quadratic function in $t$ which can be minimized.
Hint
You face an optimization problem which involves three variables $x$, $y$ and $z$ and two linear equality contraints.
I think that the simplest way is to express $y$ and $z$ as a function of $x$ using the constraints. In your case, assuming that $g_1 = x - y + 2z - 1 = 0$ and $g_2 = 2x + y - z - 1 = 0$ (which is not what is written in the top of your post : $1$ instead of $4$ as the constant term in $g_1$); this gives $y=3-5 x$ and $z=2 - 3 x$. If you replace these expressions in the objective function $$ F = D^2 = x^2 + y^2 + z^2 $$ you then have $$F=x^2+(3-5 x)^2+(2-3 x)^2=35 x^2-42 x+13$$ which is minimum when $F'=0$. Compute the corresponding value of $x$ and from there back calculate the corresponding values of $y$ and $z$.
I am sure that you can take from here.
