!****************************************************************************!
!Programme de calcul de l'exercice 1 de la la feuille 6                      !
!Comparaison des resultats donnes par                                        !
!*le schema d'Euler explicite                                                !
!*le schema de Taylor ordre 2                                                !
!*le schema de Taylor ordre 3                                                !
!*le schema de Runge-Kutta d'ordre 2                                         !
!****************************************************************************!
program comparaison
  implicit none
  real                    ::h           !le pas de temps
  real                    ::y1          !la solution exacte
  real                    ::y2,y3,y4,y5 !la valeur des differents schemas
  real                    ::t           !le temps
  integer                 ::n           !nombre de pas de temps
  integer                 ::i           !l'increment courant


  !le nombre de subdivisions
  print*,'rentrer le nombre de points'
  read*,n

  !calcul du pas de temps
  h=4.*atan(1.)/n


  !initialisation
  y1=1.
  y2=1.
  y3=1.
  y4=1.
  y5=1.
  t=0.

  open(unit=10,file='resultats.dat',status='unknown')
  write(10,*) '#temps  exacte euler  taylor2 taylor3 RK2'
  write(10,200) t,y1,y2,y3,y4,y5

  do i=1,n
     y2=y2+h*phi2(t,y2,h)
     y3=y3+h*phi3(t,y3,h)
     y4=y4+h*phi4(t,y4,h)
     y5=y5+h*phi5(t,y5,h)
     t=t+h
     y1=exp(1.-cos(t))
     write(10,200) t,y1,y2,y3,y4,y5
  enddo

  close(10)


  print*,'fin du programme'
  stop


200   format (6(E16.7,1x))


contains

  function phi2(t,y,h) result(phi)
    !schema d'Euler explicite
    implicit none
    real            ::t,y,h
    real            ::phi
    
    phi=y*sin(t)
    
  end function phi2


  function phi3(t,y,h) result(phi)
    !schema de Taylor, ordre 2
    implicit none
    real            ::t,y,h
    real            ::phi
    
    phi=y*sin(t)+0.5*h*(y*cos(t) + y*sin(t)*sin(t))
    
  end function phi3

  function phi4(t,y,h) result(phi)
    !schema de Taylor, ordre 3
    implicit none
    real            ::t,y,h
    real            ::phi
    
    phi=y*sin(t)+0.5*h*(y*cos(t) + y*sin(t)*sin(t)) & 
         &+(h**2)/6.*(-y*sin(t) +3*y*sin(t)*cos(t)+y*(sin(t))**3)
    
  end function phi4
  

  function phi5(t,y,h) result(phi)
    !schema de Runge Kutta, ordre 2
    implicit none
    real            ::t,y,h
    real            ::phi
    
    phi=sin(t+0.5*h)*(y+0.5*h*sin(t)*y)
    
  end function phi5
end program comparaison
