program rk4main

implicit none
integer :: j, n, nprint
real :: h, t, t0, tn, y, y0
real :: yrk4, k1, k2, k3, k4
character(len=10) :: output
character(len=1) :: check

problem_loop: do

   write (*,*) 'Enter initial values t0 and y0:'
   read (*,*) t0, y0
   write (*,*) 'Enter number of steps n and print frequency nprint:' 
   read (*,*) ...
   write (*,*) 'Enter final value of t:'
   read (*,*) ??
   h = ...
   t = ??
   y = ??

   ! Name and open output file:
   write (*,*) 'Enter output file name:'
   read (*,20) output
   20 format(a)
   open (unit=8, file=output, status='unknown')

   ! Header:
   write (8,30)
   30 format (//,'      t          y')
   write (8,40) t0, y0

   time_stepping_loop: do j = 1,?
      call rk4
      y = ...
      t = ...
      if (mod(j,nprint)==0) then
         write(8,40) ...
40       format(2e14.5)
      endif
   enddo time_stepping_loop

   close(8) 

   write(*,50)
50 format(//,'Do you wish to continue ? (y/n) :')
   read(*,60) check
60 format(a)
   if(check/='Y'.and.check/='y') stop
   
enddo problem_loop

stop

contains
!-----------------------------------------------

function f(t,y) result(ydash)

! Differential equation.
implicit none
real, intent(in) :: t, y
real :: ydash

ydash = ...

return
end function f
!----------------------------------------------------------------------

subroutine rk4
   ! Purpose:  To implement the classical fourth-order Runge-Kutta
   ! method for the system of M first-order differential equations
   ! y'=f(t,y) over one interval of width h.  Thus subroutine rk4
   ! approximates y(t+h) with yrk4, given y=y(t). 

k1 = h*f(t,y)
k2 = h*f(t + h/2.,y + k1/2.)
k3 = h*f(t + h/2.,y + k2/2.)
k4 = h*f(t+h,y+k3)
yrk4 = y + (k1 + 2.*k2 + 2.*k3 + k4)/6.

return
end subroutine rk4
!----------------------------------------------------------------------

end program rk4main