The most interesting part of the Euler-Lagrange equation is $\sum_i (L_p)_{x_i}$, which is the divergence of
something. Of course you know how this leads to $\Delta u$: it's the divergence of $Du$.
What field should we take so that its divergence involves $D\varphi\cdot Du$? There are two choices:
$u\,D\varphi$ and $\varphi\, Du$. Let's try them:
$$\operatorname{div}(u\,D\varphi) = D\varphi\cdot Du + u\,\Delta\varphi$$
$$\operatorname{div}(\varphi\,Du) = D\varphi\cdot Du + \varphi\,\Delta u$$
I think the second option looks more attractive, since $\Delta u$ should appear in the equation,
and $u$ by itself should not. Problem is, we have undifferentiated $\varphi$ in front of the Laplacian,
and we shouldn't.
This is where the exponential trick comes in. It's used all the time in ODE to get rid of undifferentiated functions.
For example, to solve $u'=3u$ we can write $u=e^w$ because then the equation become$e^w\,w'=3e^w$.
Here $e^w$, formerly known as $u$, cancels out, leaving a very simple equation $w'=3$.
Let's try the same idea here: instead of $\varphi$ write $e^{\varphi}$.
$$\operatorname{div}(e^\varphi\,Du) = e^\varphi D\varphi\cdot Du + e^\varphi\,\Delta u
= e^{\varphi}(D\varphi\cdot Du + \,\Delta u) $$
Good, but we need a minus sign there. So,
$$\operatorname{div}(e^{-\varphi}\,Du) = -e^{-\varphi} D\varphi\cdot Du + e^{-\varphi}\,\Delta u
= e^{-\varphi}(-D\varphi\cdot Du + \,\Delta u) $$
It remains to create $e^{-\varphi}f$ in the equation, which I leave for you to do.
OK, here it is: $L(p,z,x)=e^{-\varphi(x)}\left(\frac12|p|^2-f(x)z\right)$. Check that it works.
Best Answer
This was interesting to work on. I will humbly post my work since you don't have any other responses, even though my result has not been proven rigorously.
My conclusion is that the most general form of a Lagrangian whose Euler-Lagrange equations are generically linear-elliptic is $$L(Du,u,x)=\frac{cu^2}{2}-fu-\text{Tr}\left(A(\nabla u \otimes \nabla u) \right) $$
where $A$ is smooth and positive definite (satisfying ellipticity) and the $\otimes$ is the vector outer product to form a matrix of products of first partial derivatives. The corresponding EL equation is $$\sum_{i,j}{\bar a_{i,j}\partial_{i,j}u}+\sum_k{(\nabla\cdot \bar a_{.,k})\partial_{k}u}+cu=f$$
where I'm using $\bar a_{i,j}=A_{i,j}(1+\delta_{i,j})$ because you get a factor of 2 on the diagonal when you apply the Euler-Lagrange operator to the diagonal elements $A_{i,i}(\partial_{i}u)^2$. The divergence is of the $k$th column of the matrix $\bar a_{i,j}$.
I came to this conclusion based on the observations (again, not rigorously proved, but sound enough for me to want to post):
Any other form that is quadratic in $Du$ can be written in the form of the trace for some choice of matrix $A$. For example, $(Du)^2 = \text{Tr}(I_{n}(\nabla u \otimes \nabla u))$
No other quadratic forms (or lower) in $Du,u$ give non-trivial contributions that do not have terms in the coefficients $a_{i,j}$ that are independent of $u$. I tried them all.
Higher order derivatives and terms of degree higher than quadratic also don't give nontrivial contributions that are linear and second order without terms in the coefficients that are independent of $u$. I tried several and I think the pattern holds.
So, to answer your question directly, a completely generic linear-elliptic PDE can't be the corresponding Euler-Lagrange equation for a particular Lagrangian. The coefficients $b_j$ have to already be the divergence of the fields described by the column vectors of the matrix $\bar a_{i,j}$.
If you check the solution to the related problem you linked, you'll see that this is the case. The exponential factor $e^{-\varphi(x)}$ ensures that the desired coefficients $D\varphi$ are the divergence (derivatives) of the corresponding (diagonal) matrix representation.