PROGRAM main
IMPLICIT NONE
INTEGER,PARAMETER :: Nx=100000
INTEGER :: i
REAL :: y,x,alpha,y0
REAL :: Dx,y_old,y_new,y_ast
REAL :: rhs,x_max,x_min
REAL :: error,err_max

error = 0.
err_max = 0.

call input(x_min,x_max,alpha,y0)
Dx=(x_max-x_min)/float(Nx)
y_old=y0

call write_sol(x_min,y0,y0)

do i=1,Nx
x=x_min+Dx*i
CALL exact_sol(y,x,y0,alpha)

y_ast = y_old + Dx*rhs(y_old,alpha)
y_new = y_old + Dx*.5*(rhs(y_old,alpha)+rhs(y_ast,alpha))


error = abs(y_new-y)
err_max=max(err_max,error)

call write_sol(x,y,y_new)

y_old = y_new
enddo

WRITE(*,*) 'ERRORE MAX',Dx,err_max
WRITE(*,*) 'Fine della simulazione' 

STOP
END PROGRAM main
!_____________________________________________________________
!_____________________________________________________________
SUBROUTINE exact_sol(y,x,y0,a)
IMPLICIT NONE
REAL,INTENT(IN):: y0,x,a
REAL,INTENT(OUT):: y

y = y0*exp(-x*a)

RETURN
END SUBROUTINE exact_sol
!_____________________________________________________________
SUBROUTINE input(x_min,x_max,alpha,y0)
IMPLICIT NONE
REAL,INTENT(OUT):: y0,alpha,x_min,x_max

OPEN(UNIT=1,FILE='input.dat',STATUS='old',ACTION='read')
READ(1,*) x_min
READ(1,*) x_max
READ(1,*) y0
READ(1,*) alpha
CLOSE(1)

RETURN
END SUBROUTINE input
!_____________________________________________________________
FUNCTION rhs(y,a)
IMPLICIT NONE
REAL,INTENT(IN) :: a,y
REAL :: rhs

rhs=-a*y

RETURN
END FUNCTION rhs
!_____________________________________________________________
SUBROUTINE write_sol(x,ye,ya)
IMPLICIT NONE
REAL,INTENT(IN)::x,ye,ya

OPEN(UNIT=1,FILE='output.dat',STATUS='unknown',ACTION='write',POSITION='append')
WRITE(1,*) x,ye,ya
CLOSE(1)

RETURN
END SUBROUTINE write_sol
!____________________________________________________________