!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
            integer, parameter :: N=10                             !The number of points to use in the mesh
            double precision, dimension(0:N) :: y                  !Array used to store the solution
            double precision, parameter :: h = (b-a)/N             !This is the distance between the points
      end module RKvars

      program solveRK
            use RKvars
            double precision t
            integer count
            call RK4()
            print *,"The results of the Fourth Order Runge-Kutta Solution were:"
            do count=0,N
                  t=count*h
                  print *, "At t=", t, " the solution was:  ", y(count)
            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()
            use RKvars
            double precision, parameter :: h2 = h/2, h6 = h/6      !These are constants declared for convenience
            double precision :: K1, K2, K3, K4, th2                !The K's used in RK, and th2 declared for convenience
            double precision t
            integer count
            y(0)=init;
            do count=0, N-1                                        !Iterate through for each time step
                  t=a + (count + 1)*h
                  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) + h*K3)
                  y(count + 1) = y(count) + H6*(K1 + 2*(K2 + K3) +K4)
            end do
      end subroutine RK4