MODULE Starter 
  !Symbolic names for kind types of single- and double-precision reals: 
  INTEGER, PARAMETER :: SP = KIND(1.0) 
  INTEGER, PARAMETER :: DP = KIND(1.0D0) 
	! parameters used in program, such as number of equations in the system
	INTEGER,PARAMETER :: iter_max = 1000,no_of_equations=2
END MODULE starter

MODULE parameters

	USE Starter, ONLY : wp=>dp,no_of_equations
  IMPLICIT NONE
! declare variables used in the program
	INTEGER :: no_of_steps
  REAL(wp) :: x_start,x_finish,h,kappa
	REAL(wp) :: y_initial(no_of_equations)

CONTAINS
! read parameters from file
  SUBROUTINE input_parameters(filename)
			! declare number of integer and real parameters
      INTEGER,PARAMETER :: no_of_ints=1,no_of_reals=5

      INTEGER :: the_ints(no_of_ints),i
      real(wp) :: the_reals(no_of_reals)

      CHARACTER(LEN=*) :: filename

      OPEN(UNIT=9,file=filename)

      DO I = 1,no_of_ints
         READ(9,*)
         READ(9,*)the_ints(I)
      END DO
      DO i = 1,no_of_reals
         READ(9,*)
         READ(9,*)the_reals(i)
      ENDDO

      no_of_steps = the_ints(1);x_start=the_reals(1)
      x_finish=the_reals(2);y_initial(1)=the_reals(3)
      y_initial(2)=the_reals(4); kappa = the_reals(5)
      h = ( x_finish-x_start ) /DBLE(no_of_steps)

  END SUBROUTINE input_parameters

END MODULE parameters

MODULE funcs
	! declare modules
	USE starter, ONLY : wp=>dp,no_of_equations
  USE parameters
  IMPLICIT NONE

CONTAINS
	! function for which we solve y' = func(x,y)
	FUNCTION func(x,y)
		! declare input variables
		REAL(wp) , INTENT(IN) :: x,y(no_of_equations)
		! declare function variable
		REAL(wp) :: func(no_of_equations)
		! Here the function represents the ODE:
		!     y'' + kappa y' + x y = 0
		!     y = y_1 ; y' = y_2
		func(1) = y(2)                     !  dy_1/dx = y_2
		func(2) = - kappa * y(2) - x*y(1)  !  dy_2/dx = - kappa * y_2 - x * y_1

	END FUNCTION func

END MODULE funcs

MODULE methods
	! declare modules
	USE starter, ONLY : wp=>dp,no_of_equations
  USE funcs
  IMPLICIT NONE

CONTAINS
	! function to implement Eulers method given some function y'=func(x,y)
	FUNCTION Eulers_method(x,y,h)
		! declare input variables
		REAL(wp) , INTENT(IN) :: x,y(no_of_equations),h
		!declare output variables
		REAL(wp) :: Eulers_method(no_of_equations)
		
		Eulers_method = y + h*func(x,y)

	END FUNCTION Eulers_method

END MODULE methods

PROGRAM euler
  ! declare modules, and what we are using out of them
	USE starter, ONLY : wp=>dp,no_of_equations
  USE parameters
	USE methods, ONLY : method => Eulers_method
  IMPLICIT NONE
  ! local variables
  INTEGER :: i
  REAL(wp) :: x,y(no_of_equations)

! Input parameters from file
  CALL input_parameters("input.in")
! Open a file to write to
	OPEN(unit=10,file="euler.dat")
! Initialise x and y
	y = y_initial
	x = x_start
	! print to file
	WRITE(10,'(3(E15.8,1X))')x,y
	! run through algorithm for i = 1,n
	DO i=1,no_of_steps
		! update y = y + ...
		y = method(x,y,h)
		! update x
		x = x + h
		! print to file
		WRITE(10,'(3(E15.8,1X))')y_initial(2),x,y
	 
	END DO
	! close file
  CLOSE(unit=10)
	
	PRINT *," Value Y_1(",x_finish,") ::=",y(1)

END PROGRAM euler