You are assuming a first order autonomous dynamic $\dot x=f(x)$ and train the network $g$ so that for some fixed exact solution
\begin{align}
x(t)+g(x(t))=x(t+Δt)&=x(t)+x'(t)Δt+\frac12x''(t)Δt^2+...
\\&=x(t)+f(x(t))Δt+\frac12f^{[1]}(x(t))Δt^2+...\\
f^{[1]}(x)&=f'(x)f(x).
\end{align}
For the purpose of the task, $f$ and $g$ are unknown after construction of the training set.
Now you change the dynamic so that
$$
x_{k+1}=x_k+\frac1Pg(x_k)
$$
for sampling points with spacing $\frac1PΔt$.
If $x_0=x(t)$, then
\begin{align}
x(t+Δt)\approx x_P&=x_0+\frac1P\sum_{k=0}^{P-1}g(x_k)\\
&=x_0+\frac1P\sum_{k=0}^{P-1}f(x_k)Δt+\frac12f^{[1]}(x_k)Δt^2+...\\
&=x_0+\frac1P\sum_{k=0}^{P-1}[f(x_0)+f'(x_0)(g(x_0)+...+g(x_{k-1})]Δt+\frac12f^{[1]}(x_0)Δt^2+...\\
&=x_0+f(x_0)Δt+\frac{P-1}2f^{[1]}(x_0)Δt^2+\frac12f^{[1]}(x_0)Δt^2
\end{align}
using that $g(x)=O(Δt)$ and omitting terms of higher than quadratic degree.
As one can see, there is an excess of $\frac{P-1}2f^{[1]}(x_0)Δt^2+...$ to $x(t+Δt)$. Globally this sums up to a first-order error, not better than the Euler method. This is bad, as the assumption is that with step size $Δt$ the points sample the exact solution.
In all, you would be better served with some kind of interpolation. Compute the sequence $x_k$ with the original step size $Δt$ and use some kind of interpolation. For instance to get values inside the interval $[t_k,t_{k+1}]$,
- take the sequence $x_{k-1},x_k,x_{k+1},x_{k+1}$ and use cubic Lagrange/Newton interpolation, or
- assign derivative values via central difference quotients at $x_k,x_{k+1}$ and use, likewise cubic, Bezier/Hermite interpolation.
The first variant might be easier, the second has more consistent derivative values at the nodes. The difference should be small.
