What techniques exist for bounding the error in a solution, to an IVP, produced by RK4? The below problem is intended to contextualize the question.
Problem
The $x$, $y$ and $z$ axes of a coordiante system point to (are positive in the direction of) east, north and up, respectively. A man standing at the origin of this coordinate system throws a ball with an initial velocity of 25 m/s from a height of 1.4 m. The throw is made at a 30° angle in the -plane (vertical plane) and a 0° angle in the -plane (horizontal plane). In other words, the ball is thrown at a 30° angle (in the vertical plane) along the positive $x$-axis.
The differential equations describing the ball's trajectory are
$$ \ddot{x}=-q \dot{x}, \: \ddot{y}=-q(\dot{y}-a(z)), \: \ddot{z}=-9.81-q \dot{z}, \: \text {where} \: q=c \sqrt{\dot{x}^{2}+(\dot{y}-a(z))^{2}+\dot{z}^{2}}. $$
The air resistance coefficient $c$ is a function of the ball's radius and mass. For the ball in this problem, $c = 0.070$. The wind speed is 7 m/s on the ground, and increases with height according to $a(z) = 7 + 0.35z$.
Given the above conditions, I want to compute the ball's point of impact with an accuracy of at least $\pm10^{-4}$ m for each coordinate.
Soltuion
The above differential equations can be solved numerically by using an ODE solver such as Euler's method or RK4. I have solved them using RK4. However, to apply RK4, the equations need to be rewritten as a first-order system. This is done below.
\begin{eqnarray} s_1 = x \quad s_3 = y \quad s_5 = z \\ s_2 = \dot{x} \quad s_4 = \dot{y} \quad s_6 = \dot{z} \end{eqnarray}
\begin{equation} \begin{bmatrix} \dot{s_1} \\ \dot{s_2} \\ \dot{s_3} \\ \dot{s_4} \\ \dot{s_5} \\ \dot{s_6} \end{bmatrix} = \begin{bmatrix} s_2 \\ -qs_2 \\ s_4 \\ -q(s_4 - a(s_5)) \\ s_6 \\ -9.810 - qs_6 \\ \end{bmatrix}. \end{equation}
Now, by choosing
\begin{equation} \mathbf{s}_0 = \begin{bmatrix} 0 \\ 25\cos \pi/6 \\ 0 \\ 0 \\ 1.4 \\ 25\sin \pi/6 \end{bmatrix} \end{equation} as the starting condition and setting the step-size to $h = 0.00001$ s, the RK4 solver returns the solution $$(x, y, z) = (18.122755453098357, \: 9.016849403591824, \: -1.7008908075912395 \times 10^{-5}).$$
Question
What techniques can be used to bound the error in the computed point of impact?