! Euler's Method. A simple example

MODULE Starter 
  !Symbolic names for kind types of single- and double-precision reals: 
  INTEGER, PARAMETER :: DP = KIND(1.0D0) 
  REAL(DP),PARAMETER :: Pi = 3.141592653589793238462643383279502884197_dp 
END MODULE starter

MODULE parameters
  
  USE Starter, ONLY : wp=>dp
  IMPLICIT NONE
  ! declare parameters, not to be changed outside the routine
  INTEGER :: no_of_points_x,no_of_points_y
  REAL(wp) :: x_start,x_finish,y_start,y_finish
  
 
END MODULE parameters

MODULE funcs

  USE starter, ONLY : wp=>dp
  USE parameters
  IMPLICIT NONE

CONTAINS

  FUNCTION my_func_1(x,y)
    
    REAL(wp) , INTENT(IN) :: x,y
    REAL(wp) :: my_func_1
    ! Here the function is
    !     f(x,y) = sin(x)sin(y)
    my_func_1 = sin(x) * sin(y)
    
  END FUNCTION my_func_1

  FUNCTION my_func_2(x,y)
    
    REAL(wp) , INTENT(IN) :: x,y
    REAL(wp) :: my_func_2
    ! Here the function is
    !     f(x,y) = sin(x)cos(y)
    my_func_2 = sin(x) * cos(y)
    
  END FUNCTION my_func_2

END MODULE funcs

PROGRAM euler
  
  USE starter, ONLY : wp=>dp,pi
  USE parameters
  USE funcs
  IMPLICIT NONE
  
  INTEGER :: i,j
  REAL(wp) :: dx,dy,x,y

  ! set up initial conditions
  no_of_points_x = 100
  no_of_points_y = 100
  x_start = 0.;x_finish = 2.*pi;
  y_start = 0.;y_finish = 2.*pi;
  dx = ( x_finish-x_start ) /DBLE(no_of_points_x)
  dy = ( y_finish-y_start ) /DBLE(no_of_points_y)
  
  ! open a file to output results to
  OPEN(UNIT=10,FILE="results_twodim.dat")
  
  ! run through steps
  DO i=0,no_of_points_x
     x = dble(i)*dx
     DO j=0,no_of_points_y
	     y = dble(j)*dy	
    	 ! print values
    	 WRITE(10,'(4(e15.8,1x))')x,y,my_func_1(x,y),my_func_2(x,y)
   	 END DO
   	 WRITE(10,*)
  END DO
  
  !close file
  CLOSE(UNIT=10)
  
END PROGRAM euler