Suppose your differential equations are:
$$
\dot{x} = f_x(x,y,z)\\
\dot{y} = f_y(x,y,z)\\
\dot{z} = f_z(x,y,z),
$$
and you know $x$ and $y$ at (ordered) times $t_1, t_2, …$.
If $t_{i+1}-t_i$ is sufficiently small for some $i$ (preferably all of them), you have two ways to estimate the derivative of the first variable at $:
$$
\begin{alignat*}{1}
\dot{x}\left( \frac{t_i+t_{i+1}}{2}\right )
&\approx \frac{x(t_{i+1})-x(t_i)}{t_{i+1}-t_i} \\
\dot{x}\left( \frac{t_i+t_{i+1}}{2}\right )
&\approx f_x \left(
\frac{x(t_i)+x(t_{i+1})}{2},
\frac{y(t_i)+y(t_{i+1})}{2},
z\left( \frac{t_i+t_{i+1}}{2}\right )
\right).
\end{alignat*}
$$
If you equate the right-hand sides, the only unknown in this equation is $z\left( \frac{t_{i+1}+t_i}{2}\right )$, so all you have to do is to solve it for this term. You can then apply the same approach to $y$ to get a second estimate and refine your result.
(While I used the middle of the interval, you can apply the same approach to the beginning or end as well. This has less accuracy but may be more comfortable and even faster, as it saves one step in the following.)
Once you have this estimate for $z\left( \frac{t_{i+1}+t_i}{2}\right )$, you can integrate backwards and forwards to the known time points and refine your result. The best strategy here depends on the accuracy of your data and how your times $t_i$ are distributed. If $t_1$ is not your initial time, first obtain an accurate estimate of $z(t_1)$ and then integrate backwards to your initial time. This may not be much different from what you have been doing, but you should have a better initial guess.
Be aware of the possibility that your backwards differential equation may be unstable (e.g., for strongly dissipative systems), in which case you may want to avoid backwards integrations as far as possible in the above (or only perform them for very short time steps), in which you are essentially back at your method. Also, if $t_1$ is much later than your initial time (with respect to instability time scales), you have the more fundamental problem that the solution to your problem is overly sensitive on measurement and numerical errors – every solution to your problem would also pose a method to backwards-integrate unstable differential equations.