!HPF program for solving ordinary differential equations using Adams Fourth-Order Predictor-Corrector method

      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                        !Nused and hused pass information about current solution
            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 adams(N, h)                                 !Find the "exact" solution
            skip=N/10                                        !Only certain points are needed
            reference=y(skip:N:skip)                         !Use this solution as the reference for error calculations
            do count=1,4                                     !repeat for each time step:
                  hused=steps(count)                               !get the new time step,
                  Nused=(b-a)/hused                                !figure out how many points are needed,
                  skip=Nused/10                                    !and calculate the spacing for the output points
                  call adams(Nused, hused)                         !calculate the solution
                  error(count)=sum((reference-y(skip:Nused:skip))**2)    !total the cumulative error
            end do
            print *, "The total cumulative errors are:"             !print the results
            do count=1,4
                  print *, "Step size:  ", steps(count)
                  print *, "with a cumulative error of:  ", error(count)
                  print *, ""
            end do
      end program solveRK

      function yprime(y, t) result (value)                   !return the value of the derivative at y and t
            double precision :: y, t, value
            value=y - (t*t) + 1
      end function yprime

      subroutine RK4(hused, limit)
            use RKvars
            integer :: limit                                 !The number of points to find
            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, limit                                !Iterate through for each time step
                  t=a + (count)*hused                        !calculate the time used in all the following calculations
                  th2 = t + h2                               !calculate the Kn's
                  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)   !Update the solution
            end do
      end subroutine RK4

      subroutine adams(Nused, hused)
            use RKvars
            double precision :: hused                             !Passes the time step used
            integer :: Nused                                      !Passes the number of elements in y to use
            double precision, dimension(4:N) :: t                 !Array containing the times after the 4th step
            double precision :: h24                               !Used to eliminate repeated floating point operations
            integer count                                         !Dummy indexing variable
            t(4:Nused)=hused*[4:Nused]+a                          !set up the values for the times
            h24=hused/24
            call RK4(hused, 3)                                    !Use RK4 to get the first four points
            do count=4,Nused
                  y(count)=y(count-1) + h24*( 55*yprime(t(count-1), y(count-1)) &    !Make the initial guess
                                             -59*yprime(t(count-2), y(count-2)) &
                                             +37*yprime(t(count-3), y(count-3)) &
                                             - 9*yprime(t(count-4), y(count-4)))
                  y(count)=y(count-1) + h24*(  9*yprime(t(count), y(count)) &
                                               +19*yprime(t(count-1), y(count-1)) &  !Make the correction
                                               - 5*yprime(t(count-2), y(count-2)) &
                                               +   yprime(t(count-3), y(count-3)))
            end do
      end subroutine adams