If one sacrifices mathematical rigour, then one can get away with not using the variational bicomplex or jet bundles.
Assuming functional derivatives are defined with respect to variations with compact support, thus one can throw away all boundary terms, we have $$ \frac{\delta}{\delta \phi^a(x)}\frac{\delta}{\delta \phi^b(y)}=\frac{\delta}{\delta \phi^b(y)}\frac{\delta}{\delta \phi^a(x)}, $$
so we can try to implement the usual constructive proof of the Poincaré's lemma into a functional context. Namely, assume that $$ E_a(x)=E_a(x)[\phi] $$ is a multi-component function that depends on the field $\phi^a(x)$ in some way (usually through jets, but this approach is applicable to more general contexts), then the integrability condition for variationality is $$ \frac{\delta E_a(x)}{\delta\phi^b(x)}-\frac{\delta E_b(x)}{\delta\phi^a(x)}=0. $$ If this satisfied then locally there is a functional $S[\phi]$ such that $ E_a(x)=\frac{\delta S}{\delta\phi^a(x)}$.
The functional can be constructed explicitly by $$ S[\phi]=\int\mathrm d^nx\int_0^1\mathrm d\lambda E_a(x)[\lambda\phi]\phi^a(x). $$
The Reader may verify this by calculating the functional derivative of this expression, and using the integrability condition above.
This is just a functional analogue of the local exactness of the first degree part of the de Rham complex on finite dimensional manifolds. Namely, if $Q_\mu(x)$ are $n$ functions of $n$ variables, and it satisfies $\partial_\mu Q_\nu=\partial_\nu Q_\mu$, then at least locally (in a contractible domain) $$ Q_\mu=\partial_\mu f, $$ where $$ f(x)=\int_0^1\mathrm d\lambda Q_\mu(\lambda x)x^\mu. $$
Further remarks:
Adjustments might have to be made if the domain is contractible but does not contain 0.
The homotopy formula given in this answer agrees with that of say Anderson's The Variational Bicomplex (4.3 formula, http://deferentialgeometry.org/papers/The%20Variational%20Bicomplex.pdf)
If a Lagrangian is order $r$, then the EL equations will in general be order $2r$, although accidental degeneracies can occur which reduce the order. This procedure in general will turn an equation of motion of order $2r$ into a Lagrangian of order $2r$. Integration by parts have to be performed to bring the Lagrangian into lower order form if possible.
This constitutes a partial solution of the inverse problem to the calculus of variations. It does not however provide a general solution. If an EoM is not variational, it might be possible to use variational integrating factors, artificially increase the variable count etc. to make it variational. Moreover, if given a differential equation, it might be difficult to ascertain what the fundamental dynamical variable is (example $F_{\mu\nu}$ vs $A_\mu$ in electrodynamics).