program ode_avanzato
implicit none
integer, parameter :: n=2 ! n definisce il grado dell'equazione differenziale,
                          ! ovvero il numero di equazioni equivalenti del primo
						  ! ordine
integer            :: i, j, m_passi
real               :: h, x, x_max, x_min
real, dimension(n) :: y, yout
external           :: derivate
write(*,*) 'Soluzione dell''equazione d(2)y/dx(2) = 4*y con condizioni'
write(*,*) 'iniziali y(0)=1, y''(0)=2 (soluzione analitica y=exp(2*x))'
write(*,*)
!
! valori iniziali per l'esempio riportato in derivate
!
  x    = 0.0    ! x0
  y(1) = 1.0    ! y(x0)
  y(2) = 2.0    ! z(x0)=y'(x0)
  write(*,'(1x,t10,a,t20,a,t31,a,t47,a)') 'x','y calc','z calc=dy/dx','y analitica'
  write(*,*) x,y(1), y(2), exp(2*x)
!
!  impostazione del passo 
!
x_max = 0.1
x_min = x
m_passi = 30
h = (x_max - x_min)/m_passi   
!
! calcolo della soluzione dell'equazione differenziale
!
do i=1,m_passi
  call rk4(y,n,x,h,yout,derivate)
  x = x + h
  write(*,*) x,yout(1), yout(2), exp(2*x) 
  y = yout
enddo

end program ode_avanzato

!---------------------------------------------------------------------------
! RK4 calcola la soluzione numerica di un sistema di equazioni differenziali
! del primo ordine, utilizzando il metodo di Runge-Kutta del 4to ordine.
! x rappressenta il valore iniziale della variabile indipendente
! y rappressenta il vettore delle soluzioni, prima dell'avanzamento
! n rappressenta il numero di equazioni del primo ordine da rissolvere
! yout riporta la soluzione numerica calcolata (per tutte le equazioni)
! h rappressenta il passo
! derivs denota il nome (generico) della subroutine contenente la definizione
!        del sistema di equazioni da rissolvere (il vero nome verrà definito 
!        al momento della chiamata, e dovrà essere definito EXTERNAL nel
!        programma chiamante)
!---------------------------------------------------------------------------
subroutine rk4(y, n, x, h, yout, derivs)
implicit none
integer, intent(IN)             :: n       
real, dimension(n), intent(IN)  :: y       
real, intent(IN)                :: x, h
real, dimension(n), intent(OUT) :: yout
real                            :: h6, hh, xh
real, dimension(n)              :: k1, k2, k3, k4, yt
hh = h*0.5
h6 = h/6.
call derivs(x,y,k1,n)
xh = x + hh
yt = y + hh*k1        ! NB: le operazioni riportate coinvolgono vettori
call derivs(xh,yt,k2,n)
yt = y + hh*k2
call derivs(xh,yt,k3,n)
yt = y + h*k3
call derivs(x+h,yt,k4,n)
yout = y + h6*(k1+2.0*k2+2.0*k3+k4)
return
end subroutine rk4

!---------------------------------------------------------------------------
! derivate contiene la definizione del sistema di equazioni da rissolvere 
! x rappressenta il valore iniziale della variabile indipendente
! y rappressenta il vettore delle soluzioni, prima dell'avanzamento
! dydx rappressenta il vettore di derivate che definiscono il sistema
!      di equazioni differenziali
! n rappressenta il numero di equazioni del primo ordine da rissolvere
!---------------------------------------------------------------------------
subroutine derivate(x,y,dydx,n)
implicit none
integer, intent(in)             :: n
real, intent(in)                :: x
real, intent(in),  dimension(n) :: y
real, intent(out), dimension(n) :: dydx
! In derivate si definiscono tutte le derivate
! Esempio: risolvere l'equazione d(2)y/dx(2) = 4*y, 
!          con condizioni y(0)=1, y'(0)=2. 
!          La soluzione analitica é y=exp(2*x).
!          Mediante la sostituzione dy/dx = z,
!          l'equazione originale si trasforma nel seguente
!          sistema di equazioni:
!            dz/dx = 4*y 
!            dy/dx = z
!          con condizioni iniziali y(0)=1, z(0) = 2. 
!          Utilizzando la notazione  y(1) = y, y(2) = z,
!          si ha:
!            dydx(1) = y(2)
!            dydx(2) = 4*y(1)
 dydx(1) = y(2)
 dydx(2) = 4*y(1)
return   
end subroutine derivate