I am following the conventions of https://arxiv.org/abs/math-ph/9902027

Let $M$ be a differentiable manifold, $E \to M$ a vector bundle over $M$ with fibre $F$, $J^1(E)$ the rank-one jet bundle over $M$ and $V(M)$ the bundle of densities over $M$.

Furthermore, let $L : J^1(E) \to V(M)$ be a bundle map and define the action functional

\begin{equation} F : \Gamma(E) \to \mathbb{R}, \psi \mapsto \int\limits_M \left( L \circ j^1 \right) ( \psi ) \end{equation}

where $\Gamma(E)$ is the space of sections of $E$ and $j^1 : \Gamma(E) \to J^1(E)$ denotes the rank-one jet prolongation. How would I now derive the Euler-Lagrange equations? Assuming $L$ is nice enough, I obtain

\begin{equation} \begin{aligned} \left( D F \right)_\psi \left( \phi \right) &= \int\limits_M \left. \frac{\partial}{\partial t} \right \vert_{t = 0} \left[ \left( L \circ j^1 \right) \left( \psi + t \phi \right) \right] \\ &= \int\limits_M \left[ \left( D L \right)_{j^1 \left( \psi \right)} \circ \left( D j^1 \right)_\psi \right] \left( \phi \right) = \int\limits_M \left[ \left( D L \right)_{j^1 \left( \psi \right)} \circ j^1 \right] \left( \phi \right) \end{aligned} \end{equation}

where all $D$s can be read as Fréchet derivatives because they act on topological vector spaces. Now, I would need to perform partial integration and apply appropriate boundary conditions. One could probably introduce a positive definite inner product on $F$, translate it to the fibres and by introducing a positive definite 'reference' density $\nu$ and employing a Riesz-like representation theorem to arrive at something like

\begin{equation} \begin{aligned} \int\limits_M \left[ \left( D L \right)_{j^1 \left( \psi \right)} \circ j^1 \right] \left( \phi \right) = \int_M \langle \omega_x, \phi \rangle \, \nu(x) \end{aligned} \end{equation}

for all $\phi \in \Gamma(E)$. Here $\omega_x \in \Gamma(E)$ would be defined as

\begin{equation} \left[ \left( D L \right)_{j^1 \left( \psi \right)} \circ j^1 \right] \left( \phi \right) (x) = \langle \omega_x, \phi \rangle \, \mathrm{d} \nu(x) \end{equation}

This, however, seems rather ugly because of the arbitrariness of the reference density (I could probably live with the inner product).

The Euler-Lagrange equations should then read

\begin{equation} \omega_x = 0 \end{equation}

Of course the main issue is then still to obtain an expression for $\omega_x$ which seems difficult.