c     Programme permettant de résoudre un système différentiel
c
c     dY1/dt=F1(t,Y1,Y2,Y3)
c     dY2/dt=F2(t,Y1,Y2,Y3)
c     dY3/dt=F3(t,Y1,Y2,Y3)
c
c
      PROGRAM couple_ordre1
      implicit none
      integer*4 N,i
      real*8 h,t0,y10,y20,y30

      parameter (N=1000)          ! N=(Xmax-X0)/h ici, Xmax=10
      parameter (h=0.01d0)

      real*8 T(0:N),Y1(0:N),Y2(0:N),Y3(0:N)

      common/solution/T,Y1,Y2,Y3
      common/condinit/t0,y10,y20,y30

      t0=0.d0      !Condition initiale sur le temps
      y10=1.d0     !Condition initiale sur Y1
      y20=0.d0     !Condition initiale sur Y2
      y30=0.d0	   !Condition initiale sur Y3

      call Kutta


ccccccccccccccccccccccccccc    Affichage de la solution de l'équation différentielle par Cauchy Euler
      DO i=0,N,1
         print 120,T(i),Y1(i),Y2(i),Y3(i)              !On affiche les t(i)...y3(i) dans le format donné par 120 : x=12 entier,6décimales...
  120 format ('t=',d12.6,'   N1=',d12.6,'   N2=',d12.6,'   N3=',d12.6)   !120 est un nombre aléatoire sans importance
      ENDDO

      open(10,file='Kutta.txt',form='formatted',status='unknown')
      write(10,100)(T(i),Y1(i),Y2(i),Y3(i),i=0,n)
  100 format(d12.6,'     ',d12.6,'     ',d12.6,'     ',d12.6)
ccccccccccccccccccccccccccc


      write(*,*)' '
      write(*,*)'Tapez entrer pour fermer'
      read(*,*)


      RETURN
      END program couple_ordre1
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                    Résolution de l'équation par méthode Kutta d'ordre 4
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine Kutta
      implicit none

      integer*4 N,i
      real*8 h
      parameter(N=1000)
      parameter(h=0.01d0)

      real*8 T(0:N),Y1(0:N),Y2(0:N),Y3(0:N)
      real*8 k11,k12,k13,k14,F1
      real*8 k21,k22,k23,k24,F2
      real*8 k31,k32,k33,k34,F3
      real*8 t0,y10,y20,y30
      real*8 A1,A2,A3,A4
      real*8 B1,B2,B3,B4


      common/solution/T,Y1,Y2,Y3
      common/condinit/t0,y10,y20,y30

      T(0)=t0
      Y1(0)=y10
      Y2(0)=y20
      Y3(0)=y30

      DO i=1,N,1
         k11=h*F1(T(i-1),Y1(i-1),Y2(i-1),Y3(i-1))
         k21=h*F2(T(i-1),Y1(i-1),Y2(i-1),Y3(i-1))
         k31=h*F3(T(i-1),Y1(i-1),Y2(i-1),Y3(i-1))

         A1=T(i-1)+h/2.d0
         A2=Y1(i-1)+k11/2.d0
         A3=Y2(i-1)+k21/2.d0
         A4=Y3(i-1)+k31/2.d0

         k12=h*F1(A1,A2,A3,A4)
         k22=h*F2(A1,A2,A3,A4)
         k32=h*F3(A1,A2,A3,A4)

         B1=T(i-1)+h/2.d0
         B2=Y1(i-1)+k12/2.d0
         B3=Y2(i-1)+k22/2.d0
         B4=Y3(i-1)+k32/2.d0

         k13=h*F1(B1,B2,B3,B4)
         k23=h*F2(B1,B2,B3,B4)
         k33=h*F3(B1,B2,B3,B4)


         k14=h*F1(T(i-1)+h,Y1(i-1)+k13,Y2(i-1)+k23,Y3(i-1)+k33)
         k24=h*F2(T(i-1)+h,Y1(i-1)+k13,Y2(i-1)+k23,Y3(i-1)+k33)
         k34=h*F3(T(i-1)+h,Y1(i-1)+k13,Y2(i-1)+k23,Y3(i-1)+k33)

         T(i)=T(i-1)+h
         Y1(i)=Y1(i-1)+(1.d0/6.d0)*(k11+2.d0*k12+2.d0*k13+k14)
         Y2(i)=Y2(i-1)+(1.d0/6.d0)*(k21+2.d0*k22+2.d0*k23+k24)
         Y3(i)=Y3(i-1)+(1.d0/6.d0)*(k31+2.d0*k32+2.d0*k33+k34)
      ENDDO
      RETURN
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                    Fonction f dérivée de y.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real*8 function F1(t,y1,y2,y3)
      real*8 t,y1,y2,y3
      F1=2.d0*y1+2.d0*t     ! Expression de la fonction F1. Penser à mettre ".d0" après chaque nombre.
      RETURN
      END FUNCTION F1

      real*8 function F2(t,y1,y2,y3)
      real*8 t,y1,y2,y3
      F2=2.d0*y1+y2+4.d0*t
      RETURN
      END FUNCTION F2

      real*8 function F3(t,y1,y2,y3)
      real*8 t,y1,y2,y3
      F3=0.d0
      RETURN
      END FUNCTION F3