program decay_euler_rk4
!
! numerical integration of the differential equation for the radioactive decay using
! Euler and Runge-Kutta-4 methods
!
    implicit none
    integer,parameter:: dp=selected_real_kind(13)
    integer :: i, j, nsteps
    real(kind=dp) :: N0, N_euler, N_rk4, N_exact, lambda
    real(kind=dp) :: t, dt, tmax
    real(kind=dp) :: k1, k2, k3, k4

    ! Parameters
    N0     = 100.0_dp
    lambda = 0.20_dp
    tmax   = 10.0_dp
    
    open(unit=10,file='decay_euler_rk4.dat',status='replace')
    open(unit=1,file='decay_euler_rk4_errors.dat',status='replace')

    dt     = 0.50_dp
    
    write(10,'(A)') '# t    Euler    RK4     Exact   '

    do j= 1, 8
    
    nsteps = int(tmax/dt)

    ! Initial condition
    t = 0.0_dp
    N_euler = N0
    N_rk4   = N0
    N_exact = N0

    write(10,'(A,F10.5)') '# dt = ', dt
    write(10,'(4F15.6)') t, N_euler, N_rk4, N_exact
    
    do i = 1, nsteps

        ! Euler
        N_euler = N_euler - lambda * N_euler * dt
        
        ! RK4 step
        k1 = dt * f(N_rk4, lambda)
        k2 = dt * f(N_rk4 + 0.5_dp*k1, lambda)
        k3 = dt * f(N_rk4 + 0.5_dp*k2, lambda)
        k4 = dt * f(N_rk4 + k3, lambda)
        N_rk4 = N_rk4 + (k1 + 2.0_dp*k2 + 2.0_dp*k3 + k4)/6.0_dp

        t = t + dt

        N_exact = N0 * exp(-lambda*t)
       
        write(10,'(4F15.6)') t, N_euler, N_rk4, N_exact

    end do

    write(1,*) dt, abs(N_euler-N_exact), abs(N_rk4-N_exact)
    dt = dt / 2
    
    write(10,*)
    write(10,*)

    end do

    close(10)
    
    close (1)
    
    print *, 'Results written to decay_euler_rk4.dat and  decay_euler_rk4_errors.dat'

contains

    real(kind=dp) function f(N, lambda)
        implicit none
        real(kind=dp), intent(in) :: N, lambda

        f = -lambda * N

    end function f

end program decay_euler_rk4