program oscyllator
  real, dimension(0:2) :: trv
  external :: oscil
  integer, parameter :: neq = 2
  ! trv = (t,x,vx)
  !
  ! initial conditions
  trv(0)= 0       ! t
  trv(1) = 1      ! x(t)
  trv(2) = 0.     ! vx(t)
  dt    = 0.2;    ! time step
  t_end = 10.0;   ! final time
  !
  do
     write(*,*) trv, cos(trv(0))
     call rk4_step(neq,trv,dt,oscil)
     if (trv(0) > t_end) exit
  end do
  !
end program oscyllator

subroutine oscil(trv,der)
  real, dimension(0:2) :: trv,der
  ! calculate:  d/dt(t, x, vx) = (1, vx, ax)
  der(0) = 1        ! dt/dt   = 1 
  der(1) = trv(2)   ! dx/dt   = vx = trv(2)
  der(2) = -trv(1)  ! d2x/dt2 = ax = -x = -trv(1)
end subroutine oscil
