program main
!====================================================================
! Initial value problem for a single 1st-order ODE
! Call Methods: Euler, Predictor-Corrector, Runge-Kutta 4th order
!====================================================================
implicit none
double precision dx, dt, xi, ti, xf, tf, tmax
integer i, key
external dx

! open files (if needed)
!open (unit=6, file='result1.dat')

! initial conditions
xi = 1.00
ti = 0.00
! "time" step and max time
dt =   0.1
tmax = 2.0
! select a method: 
! key=1 Euler, key=2 predictor-corrector, key=3 RK 4th order
key = 3

!* print the header and initial conditions
write (*,*) ' 1st order single ODE '
if (key==1) write (*,*) ' Simple Euler method  ' 
if (key==2) write (*,*) ' Predictor-Corrector  '
if (key==3) write (*,*) ' Runge-Kutta 4th order'
write (*,*) '     t         x(t)   '
write (6,100) ti, xi


do while(ti <= tmax)
  tf = ti + dt

if (key==1)  call euler1(dx,ti,xi,tf,xf)
if (key==2)  call euler1m(dx,ti,xi,tf,xf)
if (key==3)  call RK4d11(dx,ti,xi,tf,xf)

  write (6,100) tf, xf
  ti = tf
  xi = xf
end do

100 format(2f10.5)
end program main

  function dx(t,x)
!----------------------------------------------
! dx(t,x)- function dx/dt in the 1st order ODE
!----------------------------------------------
implicit none
double precision dx, x, t
  dx = (-1.0)*x     ! solution x=exp(-t)
! the carrying capacity problem
!  dx = 0.1*x-0.0002*x**2 
end function dx

  subroutine RK4d11(dx,ti,xi,tf,xf)
!==================================================
! Solution of a single 1st order ODE dx/dt=f(x,t)
! Method: 4th-order Runge-Kutta method
! Alex G. February 2010
!--------------------------------------------------
! input ...
! dx(t,x)- function dx/dt (supplied by a user)
! ti - initial time
! xi  - initial position
! tf  - time for a solution
! output ...
! xf  - solution at point tf
!==================================================
implicit none
double precision dx,ti,xi,tf,xf
double precision h,k1,k2,k3,k4

h  = tf-ti

k1 = h*dx(ti,xi)
k2 = h*dx(ti+h/2.0,xi+k1/2.0)
k3 = h*dx(ti+h/2.0,xi+k2/2.0)
k4 = h*dx(ti+h,xi+k3)

xf = xi + (k1 + 2.0*(k2+k3) + k4)/6.0

end subroutine RK4d11

  subroutine euler1m(dx,ti,xi,tf,xf)
!==================================================
! Solution of a single 1st order ODE dx/dt=f(x,t)
! Method:  Modified Euler (predictor-corrector)
! Alex G. February 2010
!--------------------------------------------------
! input ...
! dx(t,x)- function dx/dt (supplied by a user)
! ti - initial time
! xi  - initial position
! tf  - time for a solution
! output ...
! xf  - solution at point tf
!==================================================
implicit none
double precision dx,xi,ti,xf,tf

xf = xi + dx(ti,xi)*(tf-ti)
xf = xi + (dx(ti,xi)+dx(tf,xf))*0.5*(tf-ti)

end subroutine euler1m

  subroutine euler1(dx,ti,xi,tf,xf)
!==================================================
! Solution of a single 1st order ODE dx/dt=f(x,t)
! Method:  Simple Euler (predictor-corrector)
! Alex G. February 2010
!--------------------------------------------------
! input ...
! dx(t,x)- function dx/dt (supplied by a user)
! ti - initial time
! xi  - initial position
! tf  - time for a solution
! output ...
! xf  - solution at point tf
!==================================================
implicit none
double precision dx,xi,ti,xf,tf

xf = xi + dx(ti,xi)*(tf-ti)

end subroutine euler1