I am attempting to approximate the Euler-Lagrange equation using a numerical solver. I have made the starting steps as follows:
Let $L = f(v, x, t)$ be provided, where $v = \frac{dx}{dt}$.
We assume that the complete numerical 3-dimensional field in $L_{(v, x, t)}$ is provided everywhere.
From the Euler-Lagrange equation $$\frac{d}{dt}(\frac{dL}{dv}) = \frac{dL}{dx}\,,$$ a spatial discretization gives $$\frac{d}{dt}(\frac{L_{(v+\Delta v, x,t)} - L_{(v,x,t)}}{\Delta v}) = \frac{L_{(v, x+\Delta x,t)} - L_{(v,x,t)}}{\Delta x}\,.$$ Then a time discretization gives $$\frac{(\frac{L_{(v+\Delta v, x,t+\Delta t)} - L_{(v,x,t+\Delta t)}}{\Delta v}) - (\frac{L_{(v+\Delta v, x,t)} - L_{(v,x,t)}}{\Delta v})}{\Delta t} = \frac{L_{(v, x+\Delta x,t)} - L_{(v,x,t)}}{\Delta x}\,.$$ Simplifying, $$\frac{L_{(v+\Delta v, x,t+\Delta t)} - L_{(v,x,t+\Delta t)} - L_{(v+\Delta v, x,t)} + L_{(v,x,t)}}{\Delta v \Delta t} = \frac{L_{(v, x+\Delta x,t)} - L_{(v,x,t)}}{\Delta x}$$ $$L_{(v+\Delta v, x,t+\Delta t)} - L_{(v,x,t+\Delta t)} - L_{(v+\Delta v, x,t)} + L_{(v,x,t)} = \frac{\Delta v \Delta t}{\Delta x}(L_{(v, x+\Delta x,t)} - L_{(v,x,t)})$$ $$L_{(v+\Delta v, x,t+\Delta t)} - L_{(v,x,t+\Delta t)} = \frac{\Delta v \Delta t}{\Delta x}L_{(v, x+\Delta x,t)} - \frac{\Delta v \Delta t}{\Delta x}L_{(v,x,t)} + L_{(v+\Delta v, x,t)} - L_{(v,x,t)} \,.$$
I don't seem to know how to proceed after this step. What is it exactly in the last line above that can be converted into a differential equation (which can then be solved using for instance the Runge-Kutta method)?
In general, that's it. You get an implicit equation for $\Delta v$, and then use $\dot x=v$ to get $\Delta x=v\,\Delta t$. Of course, one can also keep the increments for the partial derivatives independent of the actual dynamic increments of the solution. Such a separation might even make more sense if you actually explore multi-dimensional systems where you need multiple difference quotients to evaluate all partial derivatives.
A more balanced approach is to introduce impulses $p=\frac{\partial L}{\partial v}$ which give an implicit equation for $v=\dot x$ and an explicit equation $\dot p=\frac{\partial L}{\partial x}$. Note that in most mechanical systems the Lagrangian is actually separated as $L=T(v)-P(x)$ with the kinetic term $T(v)=\frac{v^TMv}2$ giving an easy relation $p=Mv$, $v=M^{-1}p$. In the next complicated stage you get a semi-mixed $T(x,v)=\frac12v^TM(x)v$ which still gives $v=M(x)^{-1}p$.