Others give useful references that discuss what is known about the answer, but no statement of the answer itself. The relevant algebraic setting is the variational bicomplex, which is discussed in the works of Anderson and others. In this setting, there are two differentials, the horizontal differential $d_H$ (representing derivatives with respect to independent variables like $t$) and the vertical differential $d_V$ (representing variational derivatives with respect to dependent variables like $q(t)$). Each of these differentials is "de Rham-like" and they anticommute with each other, which explains the cohomological flavor of the answer. A rough statement of the answer is as follows.

A Lagrangian $L$ density gives rise to a set of Euler-Lagrange equations $E_i=0$ as follows: $$ d_V L = E_i ~ d_V q^i - d_H\theta , $$ that is, the vertical 1-form $E_i ~ d_V q^i$ is vertically exact (up to a horizontally exact term $d_H \theta$). So, it is necessary for $E_i=0$ to be the Euler-Lagrange system of some Lagrangian that $E_i ~ d_V q^i$ is closed up to a horizontally exact term, namely $$ d_V(E_i~d_V q^i) = d_H \theta' (= -d_H d_V \theta) . $$ In fact, the same condition is also sufficient, up to obstructions related to the global topology of the manifold where the dependent variables $q$ take their values. This condition was formulated already classically by Helmholtz.

However, the above statement is restrictive in that it answers the question only when $E_i=0$ is already in Euler-Lagrange form. However, there are many transformations that one can apply to the system $E_i=0$ that gives an equivalent system $F_a=0$. Given only the system $F_a=0$, is it still possible to decide whether it is equivalent to some system $E_i=0$ in Euler-Lagrange form? This is the hard inverse problem (aka the multiplier problem). The only general result that I'm aware of in that direction is this.

If there exists a form $\omega$ of vertical degree 2 and horizontal degree $n-1$, where $n$ is the number of independent variables, such that it is both horizontally and vertically closed modulo the equations $F_a=0$ (namely $d_V \omega = A^a F_a$ and $d_H \omega = B^a F_a$), then there exists (again, up to global topological obstructions) a Lagrangian density $L$ whose Euler-Lagrange equations $E_i=0$ are equivalent to a subsystem of $F_a=0$.

To my knowledge, the above observation first appeared in Henneaux (AnnPhys, 1982) for ODEs and in Bridges, Hydon & Lawson (MathProcCPS, 2010) for PDEs. The calculation demonstrating this observation is given in a bit more detail on this nLab page.

It reduces the solution of the hard inverse problem to classifying all such forms $\omega$ (corresponding to the so-called characteristic cohomology of the variational bicomplex restricted to $F_a=0$ in the corresponding degree) and checking that there exists a candidate that gives rise to a Lagrangian density whose Euler-Lagrange system $E_i=0$ is equivalent to the full system $F_a=0$. The calculation of the corresponding characteristic cohomology of the system $F_a=0$ is still non-trivial, but there exist ways of attacking it, which include Vinogradov's $\mathcal{C}$-spectral sequence mentioned in other responses.

Others give useful references that discuss what is known about the answer, but no statement of the answer itself. The relevant algebraic setting is the variational bicomplex, which is discussed in the works of Anderson and others. In this setting, there are two differentials, the horizontal differential $d_H$ (representing derivatives with respect to independent variables like $t$) and the vertical differential $d_V$ (representing variational derivatives with respect to dependent variables like $q(t)$). Each of these differentials is "de Rham-like" and they anticommute with each other, which explains the cohomological flavor of the answer. A rough statement of the answer is as follows.

A Lagrangian $L$ density gives rise to a set of Euler-Lagrange equations $E_i=0$ as follows: $$ d_V L = E_i ~ d_V q^i - d_H\theta , $$ that is, the vertical 1-form $E_i ~ d_V q^i$ is vertically exact (up to a horizontally exact term $d_H \theta$). So, it is necessary for $E_i=0$ to be the Euler-Lagrange system of some Lagrangian that $E_i ~ d_V q^i$ is closed up to a horizontally exact term, namely $$ d_V(E_i~d_V q^i) = d_H \theta' (= -d_H d_V \theta) . $$ In fact, the same condition is also sufficient, up to obstructions related to the global topology of the manifold where the dependent variables $q$ take their values. This condition was formulated already classically by Helmholtz.

However, the above statement is restrictive in that it answers the question only when $E_i=0$ is already in Euler-Lagrange form. However, there are many transformations that one can apply to the system $E_i=0$ that gives an equivalent system $F_a=0$. Given only the system $F_a=0$, is it still possible to decide whether it is equivalent to some system $E_i=0$ in Euler-Lagrange form? This is the hard inverse problem (aka the multiplier problem). The only general result that I'm aware of in that direction is this.

If there exists a form $\omega$ of vertical degree 2 and horizontal degree $n-1$, where $n$ is the number of independent variables, such that it is both horizontally and vertically closed modulo the equations $F_a=0$ (namely $d_V \omega = A^a F_a$ and $d_H \omega = B^a F_a$), then there exists (again, up to global topological obstructions) a Lagrangian density $L$ whose Euler-Lagrange equations $E_i=0$ are equivalent to a subsystem of $F_a=0$.

To my knowledge, the above observation first appeared in Henneaux (AnnPhys, 1982) for ODEs and in Bridges, Hydon & Lawson (MathProcCPS, 2010) for PDEs. The calculation demonstrating this observation is given in a bit more detail on this nLab page. (Edit: At risk of shameless self-promotion, I'll also note that I collected these observations in a self-contained paper (arXiv; JMP, 2013).)

It reduces the solution of the hard inverse problem to classifying all such forms $\omega$ (corresponding to the so-called characteristic cohomology of the variational bicomplex restricted to $F_a=0$ in the corresponding degree) and checking that there exists a candidate that gives rise to a Lagrangian density whose Euler-Lagrange system $E_i=0$ is equivalent to the full system $F_a=0$. The calculation of the corresponding characteristic cohomology of the system $F_a=0$ is still non-trivial, but there exist ways of attacking it, which include Vinogradov's $\mathcal{C}$-spectral sequence mentioned in other responses.

Others give useful references that discuss what is known about the answer, but no statement of the answer itself. The relevant algebraic setting is the variational bicomplex, which is discussed in the works of Anderson and others. In this setting, there are two differentials, the horizontal differential $d_H$ (representing derivatives with respect to independent variables like $t$) and the vertical differential $d_V$ (representing variational derivatives with respect to dependent variables like $q(t)$). Each of these differentials is "de Rham-like" and they anticommute with each other, which explains the cohomological flavor of the answer. A rough statement of the answer is as follows.

A Lagrangian $L$ density gives rise to a set of Euler-Lagrange equations $E_i=0$ as follows: \[ d_V L = E_i ~ d_V q^i + d_H\theta , \] that is, the vertical 1-form $E_i ~ d_V q^i$ is vertically exact (up to a horizontally exact term $d_H \theta$). So, it is necessary for $E_i=0$ to be the Euler-Lagrange system of some Lagrangian that $E_i ~ d_V q^i$ is closed up to a horizontally exact term, namely \[ d_V(E_i~d_V q^i) = d_H \theta' (= d_H d_V \theta) .\] In fact, the same condition is also sufficient, up to obstructions related to the global topology of the manifold where the dependent variables $q$ take their values. This condition was formulated already classically by Helmholtz.

However, the above statement is restrictive in that it answers the question only when $E_i=0$ is already in Euler-Lagrange form. However, there are many transformations that one can apply to the system $E_i=0$ that gives an equivalent system $F_a=0$. Given only the system $F_a=0$, is it still possible to decide whether it is equivalent to some system $E_i=0$ in Euler-Lagrange form? This is the hard inverse problem (aka the multiplier problem). The only general result that I'm aware of in that direction is this.

If there exists a form $\omega$ of vertical degree 2 and horizontal degree $n-1$, where $n$ is the number of independent variables, such that it is both horizontally and vertically closed modulo the equations $F_a=0$ (namely $d_V \omega = A^a F_a$ and $d_H \omega = B^a F_a$), then there exists (again, up to global topological obstructions) a Lagrangian density $L$ whose Euler-Lagrange equations $E_i=0$ are equivalent to a subsystem of $F_a=0$.

To my knowledge, the above observation first appeared in Henneaux (AnnPhys, 1982) for ODEs and in Bridges, Hydon & Lawson (MathProcCPS, 2010) for PDEs. The calculation demonstrating this observation is given in a bit more detail on this nLab page.

It reduces the solution of the hard inverse problem to classifying all such forms $\omega$ (corresponding to the so-called characteristic cohomology of the variational bicomplex restricted to $F_a=0$ in the corresponding degree) and checking that there exists a candidate that gives rise to a Lagrangian density whose Euler-Lagrange system $E_i=0$ is equivalent to the full system $F_a=0$. The calculation of the corresponding characteristic cohomology of the system $F_a=0$ is still non-trivial, but there exist ways of attacking it, which include Vinogradov's $\mathcal{C}$-spectral sequence mentioned in other responses.

Others give useful references that discuss what is known about the answer, but no statement of the answer itself. The relevant algebraic setting is the variational bicomplex, which is discussed in the works of Anderson and others. In this setting, there are two differentials, the horizontal differential $d_H$ (representing derivatives with respect to independent variables like $t$) and the vertical differential $d_V$ (representing variational derivatives with respect to dependent variables like $q(t)$). Each of these differentials is "de Rham-like" and they anticommute with each other, which explains the cohomological flavor of the answer. A rough statement of the answer is as follows.

A Lagrangian $L$ density gives rise to a set of Euler-Lagrange equations $E_i=0$ as follows: $$ d_V L = E_i ~ d_V q^i - d_H\theta , $$ that is, the vertical 1-form $E_i ~ d_V q^i$ is vertically exact (up to a horizontally exact term $d_H \theta$). So, it is necessary for $E_i=0$ to be the Euler-Lagrange system of some Lagrangian that $E_i ~ d_V q^i$ is closed up to a horizontally exact term, namely $$ d_V(E_i~d_V q^i) = d_H \theta' (= -d_H d_V \theta) . $$ In fact, the same condition is also sufficient, up to obstructions related to the global topology of the manifold where the dependent variables $q$ take their values. This condition was formulated already classically by Helmholtz.

However, the above statement is restrictive in that it answers the question only when $E_i=0$ is already in Euler-Lagrange form. However, there are many transformations that one can apply to the system $E_i=0$ that gives an equivalent system $F_a=0$. Given only the system $F_a=0$, is it still possible to decide whether it is equivalent to some system $E_i=0$ in Euler-Lagrange form? This is the hard inverse problem (aka the multiplier problem). The only general result that I'm aware of in that direction is this.

If there exists a form $\omega$ of vertical degree 2 and horizontal degree $n-1$, where $n$ is the number of independent variables, such that it is both horizontally and vertically closed modulo the equations $F_a=0$ (namely $d_V \omega = A^a F_a$ and $d_H \omega = B^a F_a$), then there exists (again, up to global topological obstructions) a Lagrangian density $L$ whose Euler-Lagrange equations $E_i=0$ are equivalent to a subsystem of $F_a=0$.

To my knowledge, the above observation first appeared in Henneaux (AnnPhys, 1982) for ODEs and in Bridges, Hydon & Lawson (MathProcCPS, 2010) for PDEs. The calculation demonstrating this observation is given in a bit more detail on this nLab page.

It reduces the solution of the hard inverse problem to classifying all such forms $\omega$ (corresponding to the so-called characteristic cohomology of the variational bicomplex restricted to $F_a=0$ in the corresponding degree) and checking that there exists a candidate that gives rise to a Lagrangian density whose Euler-Lagrange system $E_i=0$ is equivalent to the full system $F_a=0$. The calculation of the corresponding characteristic cohomology of the system $F_a=0$ is still non-trivial, but there exist ways of attacking it, which include Vinogradov's $\mathcal{C}$-spectral sequence mentioned in other responses.