Clear[EulerLagrange];
Options[EulerLagrange] = {eXpand -> False};
EulerLagrange[density_, depend_List, independ_List, options__] :=
Block[{f0, fh, \[Epsilon], w, y, x$m, expand, euler = {}, wtable},
{expand} = {eXpand} /. {options};
wtable = Table[w[i], {i, 1, Length[depend]}];
f0 = Function[x$m, y + \[Epsilon] w];
ruleg[i_] :=
b_ Derivative[n__][wtable[[i]]][independ] :> (-1)^
Plus @@ {n} HoldForm[\!\(
\*SubscriptBox[\(\[PartialD]\), \(Delete[Thread[{independ, {n}}],
0]\)]b\)];
Do[
fh = density /. depend[[j]] -> f0 /. {x$m -> independ,
y -> depend[[j]] @@ independ, w -> wtable[[j]] @@ independ};
fh = Expand[\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Epsilon]\)]fh\) /. \[Epsilon] ->
0];
fh = fh /. ruleg[j] /. wtable[[j]] @@ independ -> 1;
AppendTo[euler, fh], {j, 1, Length[depend]}];
If[Not[expand], euler = ReleaseHold[euler], euler]]
f = (\!\(
\*SubscriptBox[\(\[PartialD]\), \(x\)]\(u[x, y, z]\)\))^2 + (\!\(
\*SubscriptBox[\(\[PartialD]\), \(y\)]\(u[x, y, z]\)\))^2 + (\!\(
\*SubscriptBox[\(\[PartialD]\), \(z\)]\(u[x, y, z]\)\))^2
EulerLagrange[f, {u}, {x, y, z}, eXpand -> False]

顺便一提，VariationalMethods包有一个列拉格朗日方程的函数。