t}{2}, h_n + dh_1\frac{\Delta t}{2}, V_n + dv_1\frac{\Delta t}{2}, \nonumber\\
  && \quad q(V_n + dv_1\frac{\Delta t}{2},h_n + dh_1\frac{\Delta t}{2})) \nonumber\\
  dv_2 &=& f_2(t_n + \frac{\Delta t}{2}, V_n + dv_1\frac{\Delta t}{2},h_n + dh_1\frac{\Delta t}{2}, dh_2) \nonumber \\
  dh_3 &=& f_1(t_n + \frac{\Delta t}{2}, h_n + dh_2\frac{\Delta t}{2}, V_n + dv_2\frac{\Delta t}{2}, \nonumber\\
  && \quad q(V_n + dv_2\frac{\Delta t}{2},h_n + dh_2\frac{\Delta t}{2})) \nonumber\\
  dv_3 &=& f_2(t_n + \frac{\Delta t}{2},V_n + dv_2\frac{\Delta t}{2},h_n + dh_2\frac{\Delta t}{2}, dh_3) \nonumber \\
  dh_4 &=& f_1(t_n + \Delta t, h_n + dh_3 \Delta t, V_n + dv_3 \Delta t, \nonumber\\
  && \quad q(V_n + dv_3 \Delta t,h_n + dh_3 \Delta t)) \nonumber\\
  dv_4 &=& f_2(t_n + \Delta t, V_n + dv_3 \Delta t,h_n + dh_3 \Delta t, dh_4) \nonumber \\  
  h_{n+1} &=& h_n + \frac{\Delta t}{6} \left( dh_1 + 2dh_2 + 2dh_3 + dh_4 \right) \nonumber \\  
  V_{n+1} &=& V_n + \frac{\Delta t}{6} \left( dv_1 + 2dv_2 + 2dv_3 + dv_4 \right) \nonumber \\
  t_{n+1} &=& t_n + \Delta t \nonumber\end{aligned}$$ where the initial conditions $t_0$, $q_0$, $h(t_0,q_0)=h_0 $, $V(t_0,h_0)=V_0$ are given. The choice of the Runge-kutta method is justified by the jumps in the function $f$. As a matter of fact, this method will minimize the 