program spring_midpoint_fnf
  IMPLICIT NONE

  integer, parameter :: DP = SELECTED_REAL_KIND(14)

  integer, parameter :: nvars = 2

  real( kind=DP ), dimension(nvars) :: y, ymid, f
  real( kind=DP ) :: vinit = 5.0_DP, xinit = 0.0_DP, time=0.0_DP
  real( kind=DP ) :: en, step
  character(20) :: nchar, stepchar !strings for reading commandline
  integer :: n,i

!  Get n and step from command line
!  make sure we have both arguments
  if(command_argument_count() .ne. 2)then
    write(*,*)'Usage: spring-midpoint-fnf n step'
    stop
  endif
  call get_command_argument(1,nchar) !read first argument
  call get_command_argument(2,stepchar) !read 2nd argument

  read(nchar,*)n  ! convert string to whatever type n is (integer)
  read(stepchar,*)step 

  i = 0 ! The loop counter

  y(1) = xinit ! Initial conditions
  y(2) = vinit

  en = energy(y)
!      see function definition below

  write(*,*) time, y(1), y(2), en
  do while( i < n )
    call deriv(y,f)  !see subroutine below 
    ymid(1) = y(1) + 0.5_DP*step*f(1)
    ymid(2)= y(2) + 0.5_DP*step*f(2)

    call deriv(ymid,f)
    y(1) = y(1)+step*f(1)
    y(2) = y(2)+step*f(2)

    en = energy(y)
    time = time+step
    i=i+1

    write(*,*) time, y(1), y(2), en
  enddo

  call flush()
! On my system, for some reason, not all of the output
! was being captured when I was redirecting standard
! output to a full.  Calling this subroutine makes sure
! that anything being written out gets flushed out of
! any buffer that is holding the information temporarily.

  contains

  function energy(y)

    integer, parameter :: DP = SELECTED_REAL_KIND(14)
    real( kind=DP ) :: energy
    real( kind=DP ), dimension(*), intent(in) :: y
! dimension(*) means array will take the size of
! array the function was called with

!   recall that m=1 and k=1, y(1) is x, and y(2) is vx
    energy = 0.5_DP*(y(1)**2 + y(2)**2)
    return
  end function energy

end program

subroutine deriv(y,f)
  IMPLICIT NONE

  integer, parameter :: DP = SELECTED_REAL_KIND(14)
  real( kind=DP ), dimension(*), intent(in) :: y
  real( kind=DP ), dimension(*), intent(inout) :: f

  f(1) = y(2)
  f(2) = -y(1)

  return

end subroutine