PROGRAM main
IMPLICIT NONE
INTEGER,PARAMETER :: Nx=5000
INTEGER,PARAMETER :: nrk=4
INTEGER :: i,l
REAL :: y,x,alpha,y0
REAL :: Dx,y_old,y_new,y_ast
REAL :: rhs,x_max,x_min
REAL :: error,err_max
REAL,DIMENSION(1:nrk) :: crk1,crk2
REAL,DIMENSION(0:nrk) :: k
REAL :: rhsrk

crk1= (/0., .5, .5, 1./)
crk2= (/1., 2., 2., 1./)

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)

rhsrk=0.
y_new = y_old
k(0)=0.
do l=1,nrk
y_old = y_new + crk1(l) * Dx * k(l-1)
k(l) = rhs(y_old,alpha)
rhsrk = rhsrk + k(l)*crk2(l)
enddo

y_new = y_new + Dx/6.*rhsrk


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
!____________________________________________________________