PROGRAM Eulers_Method
!-----------------------------------------------------------------------
! Program that uses Euler's method to obtain an approximate solution 
! to a first order differential equation of the form:
!                           Y' = F(X, Y)
! Variables used are:
!   X              : current X-value 
!   X_Next         : next X-value (X + DeltaX)
!   Y              : approximate Y-value corresponding to X
!   DeltaX         : X-increment used
!   NumIterations  : number of iterations 
!   N              : counts iterations
!   F              : function in the diff. equation (internal function)
!
! Input:  Initial values for X and Y, DeltaX, and NumIterations 
! Output: A sequence of points (X, Y) that approximate the solution
!         curve 
!-----------------------------------------------------------------------

  IMPLICIT NONE
  REAL :: X, Y, X_Next, DeltaX
  INTEGER :: N, NumIterations

  ! Get given information; print table headings and initial values

  PRINT *, "Enter X0 and Y0, X-increment to use, and"
  PRINT *, "the number of values to calculate:"
  READ *, X, Y,  DeltaX, NumIterations
  PRINT *
  PRINT *, "      X          Y"
  PRINT *, "===================="
  PRINT '(1X, 2F10.5)',  X, Y 


  ! Iterate with Euler's method

  DO N = 1, NumIterations
     X_Next = X + DeltaX
     Y = Y + F(X,Y) * DeltaX
     X = X_Next
     PRINT '(1X, 2F10.5)', X, Y 
  END DO

CONTAINS

  !- F(X, Y) --------------------------------------------
  ! The function F in differential equation Y" = F(X, Y)
  !------------------------------------------------------
  
  FUNCTION F(X, Y)
  
    REAL :: F
    REAL, INTENT(IN) :: X, Y
  
    F = 2.0 * X * Y
  
  END FUNCTION F
  
END PROGRAM Eulers_Method