!
! =====> Program - P123.F90
!
        PROGRAM DifDmo
        IMPLICIT NONE
!  Program to demonstrate numerical solutions
!  to ordinary, first order differential equations.

        INTEGER :: I, N
        REAL :: X, Xf, h,        &
                Y0, Y_Eul, Y_RK4

        PRINT *, 'This is Program >> P123 = Differential Equations'
!
!
        PRINT *, 'First order ODEs   F = -y*y  for example'
        PRINT *, 'Enter final (high) x value: '
        READ *, Xf
        Print *, Xf

        PRINT *, 'Number of steps to take: '
        READ *, N
        Print *, N

!       Step size, initial X, Y (given)
        h = Xf/N
        X = 0.0
        Y0 = 1

        PRINT *, 'Step size: ', h
        PRINT 101, 'Step', 'x value', 'Euler', 'RK-4'

!       Compare Euler and Runge-Kutta 4'th order
!       at each step.
        Y_Eul = Y0
        Y_RK4 = Y0
L1:     DO I = 1, N
          Y_Eul = Euler(X, Y_Eul, h)
          Y_RK4 = RngKt4(X, Y_RK4, h)

          X = X + h

          PRINT 103, I, X, Y_Eul, Y_RK4
        END DO L1

        STOP
  101   FORMAT(//' ', A5, 6A12)
  103   FORMAT(' ', I5, 3F12.8)

        END PROGRAM DifDmo


        FUNCTION F(x, y)
!       The function to be integrated.
!       I.e. y' = dy/dx = F(x, y)

        F = x     ! Just to keep FTN90 quiet
        F = -y*y
        RETURN
        END FUNCTION F


! ALGORITHMS

        FUNCTION Euler(x, y, h)
!       Uses F(x, y) to advance the solution
!       from x to x+h.

        REAL :: x, y, h

        Euler = y + h*F(x, y)

        RETURN
        END FUNCTION Euler



        FUNCTION RngKt4(x, y, h)
!       Uses F(x, y) to advance the solution
!       from x to x+h.

        REAL ::         x, y, h,  &
                        h_hf,     &
                        x_hf,     &
                        y_hf_1,   &
                        y_hf_2,   &
                        y_1

        h_hf  = h/2
        x_hf  = x + h_hf

!       Form various estimates of the solution
!       at x+h/2 and x+h, then combine them
!       using the formula below.
        y_hf_1  = Euler(x, y, h_hf)
        y_hf_2  = Euler(x_hf, y_hf_1, h_hf)
        y_1     = Euler(x_hf, y_hf_2, h)

        RngKt4 = y + h/3 *                  &
          ((f(x, y) + f(x+h, y_1))/2 +      &
          f(x_hf, y_hf_1) + f(x_hf, y_hf_2))

      RETURN
      END FUNCTION RngKt4