The Discussion

This is a very simple, almost a self-explanatory code. I'll skip most of the variable definitions and the assert_eq bit, because it's much the same as before.

xs

is the starting point, x₀

htot

is the total step $\Delta x$ to be taken, but it's going to be split into multiple shorter steps

nstep

is the number of those shorter steps to be used

We begin by evaluating the length of a mid-point step by dividing our input value of $\Delta x$ by nstep:

  h = htot / nstep

Thus $h = \Delta{x} / n$. Then we set y_m, to y and y_n to y + dydx h and set x = x₀ + h:

  ym = y
  yn = y + h * dydx
  x = xs + h

So that at this stage y_m = y₀ and y_n = y₁, whereas x = x₁.

The next step is to evaluate the right hand side of the differential equation at x₁ and y₁:

  CALL derivs(x, yn, yout)
  h2 = 2.0_sp * h

Now we are ready to start the leapfrog. h₂ becomes 2 h and we enter the DO loop within which:

1.

we swap the roles of y_m and y_n, so that the guy that has become an endpoint as the result of the last step, becomes the midpoint point now, whereas the guy that was a midpoint before becomes a starting point now. Time, i.e., x, in the meantime becomes advanced by half-step, i.e., h:

     CALL swap(ym, yn)
     yn = yn + h2 * yout
     x = x + h

2.

having gone through these preparations we evaluate the right hand side of the equation at the new end-point that is going to become a mid-point for the next jump:

     CALL derivs(x, yn, yout)

The last step, just before leaving the subroutine, is to wind down the process by taking the average of the implicitly evaluated last value of y, i.e., ym + h * yout and explicitly evaluated value, i.e., yn, and return that average in yout:

  yout = 0.5_sp * (ym + yn + h * yout)