C		Subroutine de resolution d'une equation differentielle
C		par la methode de Rung-Kutta d'ordre 4
C		Dominique LEFEBVRE mars 2006

	  SUBROUTINE ResolutionRK4(x,y,dx,NbEquation,Derivee)
C		x    	   : abscisse
C		y()        : valeurs retournees par le systeme dif.
C		dx         : pas de calcul
C		NbEquation : nombre d'equations dans le systeme dif.
C		Derivee    : nom de la fonction decrivant le systeme dif.
	  
      IMPLICIT NONE			
	  EXTERNAL Derivee		! fonction de calcul de la derivee
C
C		Typage de l'interface de ResolutionRK4
C		
	  INTEGER NbEquation
	  REAL x, y(NbEquation), dx
C		
C		Declaration des variables locales
	  REAL pred1(NbEquation),pred2(NbEquation),pred3(NbEquation),
	 1		    pred4(NbEquation), ytemp(NbEquation), halfdx
	  INTEGER i

	  halfdx = dx/2

C		Premiere prediction
	  CALL Derivee(x,y,pred1,NbEquation)
	  DO i=1, NbEquation
		ytemp(i) = y(i) + pred1(i)*halfdx
	  ENDDO 		
C 		Seconde prediction
	  CALL Derivee(x+halfdx,ytemp,pred2,NbEquation)
	  DO i=1, NbEquation
		ytemp(i) = y(i) + pred2(i)*halfdx
	  ENDDO 
C		Troisieme prediction
	  CALL Derivee(x+halfdx,ytemp,pred3,NbEquation)
	  DO i=1, NbEquation
		ytemp(i) = y(i) + pred3(i)*dx
	  ENDDO 				
C		Quatrieme prediction		
      CALL Derivee(x+dx,ytemp,pred4,NbEquation)

C		Estimation de y pour le pas suivant
	  DO i=1, NbEquation
		y(i)=y(i)+(dx/6.0)*(pred1(i)+2.0*pred2(i)+2.0*pred3(i) + pred4(i))	  
	  ENDDO

	  END