!HPF program for solving ordinary differential equations using fourth order Runge-Kutte

      module RKvars
            double precision, parameter :: a=0.0, b=2.0              !These are the endpoints for the solution
            double precision, parameter :: init=0.5                  !This is the initial value
            double precision, parameter :: h=.001                    !The step size (in time)
            integer, parameter :: N=(b-a)/h                          !The number of points to use in the mesh
            double precision, dimension(0:N) :: y                    !Array used to store the solutions
      end module RKvars

      program solveRK
            use RKvars
            implicit none
            integer :: skip                                          !Included for convienence
            integer :: count, status, Nused                          !dummy variables, status verifies allocation of arrays
            double precision :: hused
            double precision, dimension(1:4) :: error, steps         !Arrays where error and time step, values will be stored
            double precision, dimension(1:10)  :: reference          !Stores the reference points (assumed to be exact)
            steps=[0.005, 0.010, 0.020, 0.040]                       !The step sizes required in the hw problem
            call RK4(N, h)
            skip=N/10
            reference=y(skip:N:skip)
            do count=1,4
                  hused=steps(count)
                  Nused=(b-a)/hused
                  skip=Nused/10
                  call RK4(Nused, hused)
                  error(count)=sum((reference-y(skip:Nused:skip))**2)
            end do
            print *, "The total cumulative errors are:"
            do count=1,4
                  print *, "step size:  ", steps(count)
                  print *, "with cumulative Error:  ", error(count)
                  print *,""
            end do
      end program solveRK

      function yprime(y, t) result (value)
            double precision :: y, t, value
            value=y - (t*t) + 1
      end function yprime

      subroutine RK4(Nused, hused)
            use RKvars
            integer :: Nused                                         !The number of points actually used
            double precision :: hused                                !The step size actually used
            double precision :: h2, h6                               !These are declared for convenience
            double precision :: K1, K2, K3, K4, th2                  !The K's used in RK, and th2 declared for convenience
            double precision :: t                                    !time
            integer count
            h2=hused/2.0
            h6=hused/6.0
            y=0.0
            y(0)=init;
            do count=0, Nused-1                                       !Iterate through for each time step
                  t=a + (count)*hused
                  th2 = t + h2
                  K1 = yprime(y(count), t)
                  K2 = yprime(th2, y(count) + h2*K1)
                  K3 = yprime(th2, y(count) + h2*K2)
                  K4 = yprime(t + h, y(count) + hused*K3)
                  y(count + 1) = y(count) + H6*(K1 + 2*(K2 + K3) +K4)
            end do
      end subroutine RK4