!********************************************************************
!*             Differential equations of order 1                    *
!*             by Runge-Kutta method of order 4                     *
!*                                                                  *
!*                               F90 version by J-P Moreau, Paris   *
!*                                      (www.jpmoreau.fr)           *
!* ---------------------------------------------------------------- *
!* Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc  *
!*             DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991"    *
!*             [BIBLI 03].                                          *
!* ---------------------------------------------------------------- *
!* SAMPLE RUN:                                                      *
!*                                                                  *
!* Example: integrate y'=4x*(y+sqrt(y))/1+x^2 from x=0 to x=1       *
!*                                                                  *
!*        DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER 1          *
!*                     of type y' = f(x,y)                          *
!*                                                                  *
!*   begin value x   : 0                                            *
!*   end value x     : 1                                            *
!*   y value at x0   : 1                                            *
!*   number of points: 11                                           *
!*   finesse         : 10                                           *
!*                                                                  *
!*       X            Y                                             *
!* -------------------------                                        *
!*   0.000000     1.000000                                          *
!*   0.100000     1.040400                                          *
!*   0.200000     1.166400                                          *
!*   0.300000     1.392400                                          *
!*   0.400000     1.742400                                          *
!*   0.500000     2.250000                                          *
!*   0.600000     2.958400                                          *
!*   0.700000     3.920400                                          *
!*   0.800000     5.198400                                          *
!*   0.900000     6.864400                                          *
!*   1.000000     9.000000                                          *
!*                                                                  *
!********************************************************************
PROGRAM TEST_EQDIF1

integer fi,m
real*8 xi,xf,yi

  print *,' '
  print *,' DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER 1'
  print *,'              of type y'' = f(x,y)'
  print *,' '
  write(*,"('    begin value x   : ')",advance='no')
  read *, xi
  write(*,"('    end value x     : ')",advance='no')
  read *, xf
  write(*,"('    y value at x0   : ')",advance='no')
  read *, yi
  write(*,"('    number of points: ')",advance='no')
  read *, m
  write(*,"('    finesse         : ')",advance='no')
  read *, fi

  call Equadif1(xi,xf,yi,m,fi)

end

real*8 Function fp(x, y)
  !Example: y'=4x(y+rac(y))/(1+x²)
  real*8 x,y
  if (y<0.d0) y=0.d0
  fp = 4.d0*x*(y+dsqrt(y))/(1.d0+x*x)
end

Subroutine Affiche(t,m,xi,xf)
  parameter(SIZE=50) 
  real*8 fp,t(SIZE),xi,xf,h,x
  integer i,m
  h=(xf-xi)/(m-1)
  x=xi-h
  print *,' '
  print *,'      X           Y     '
  print *,'------------------------'
  do i=1, m
    x=x+h
    write(*,10) x, t(i)
  end do
  print *,' '
  return
10 format(' ',f9.6,'   ',f9.6)
end

!***************************************************************************
!*        SOLVING DIFFERENTIAL EQUATIONS WITH 1 VARIABLE OF ORDER 1        *
!*                       of type y' = f(x,y)                               *
!* ----------------------------------------------------------------------- *
!*  INPUTS:                                                                *
!*    fp        Given equation to integrate (see test program)             *
!*    xi, xf    Begin, end values of variable x                            *
!*    yi        Begin Value of y at x=xi                                   *
!*    m         Number of points to calculate                              *
!*    fi        finesse (number of intermediate points)                    *
!*                                                                         *
!*  OUTPUTS:                                                               *
!*    t       : real vector storing m results for function y               *
!***************************************************************************
Subroutine Equadif1(xi,xf,yi,m,fi)
  parameter(SIZE=50)
  integer fi,i,j,m,ni
  real*8 t(SIZE),xi,xf,yi,a,b,c,d,h,x,y;
  if (fi < 1) return
  h  = (xf - xi) / fi / (m-1)
  y  = yi
  t(1) = yi

  do i=1, m
    ni = (i - 1) * fi - 1
    do j=1, fi
      x = xi + h * (ni + j)
      a = h * fp(x,y)
      b = h * fp(x+h/2.d0,y+a/2.d0)
      c = h * fp(x+h/2.d0,y+b/2.d0)
      x = x + h
      d = h * fp(x,y+c)
      y = y + (a + b + b + c + c + d) / 6.d0
    end do
    t(i+1) = y
  end do
  call Affiche(t,m,xi,xf)
  return
end

! End of file teqdif1.f90