!***********************************************************
!*        SOLVING DIFFERENTIAL EQUATIONS OF ORDER 1        *
!*             (Prediction-correction method)              *
!* ------------------------------------------------------- *
!* Reference:                                              * 
!*   "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc  *
!*    DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991"    *
!*    [BIBLI 03].                                          *  
!*                                                         *
!*                      F90 version by J-P Moreau, Paris   *
!*                             (www.jpmoreau.fr)           *
!* ------------------------------------------------------- *
!* SAMPLE RUN:                                             *
!* (Integrate y'=(4xy+4x*sqrt(y))/(1+x²) from x=0 to x=10  *
!*  with initial condition y(0)=1)                         *
!*                                                         *
!*  Input x begin: 0                                       *                 
!*  Input x end  : 10                                      *
!*  Input y begin: 1                                       * 
!*  Input number of points: 10                             *
!*  Input finesse: 10                                      *
!*                                                         *
!*         X               Y                               *
!*  -----------------------------                          *
!*     1.000000         8.999984                           *
!*     2.000000        80.999881                           *
!*     3.000000       360.999496                           *
!*     4.000000      1088.998512                           *
!*     5.000000      2600.996484                           *
!*     6.000000      5328.992838                           *
!*     7.000000      9800.986875                           *
!*     8.000000     16640.977767                           *
!*     9.000000     26568.964559                           *
!*    10.000000     40400.946171                           *
!*  -----------------------------                          *
!*                                                         *
!***********************************************************
PROGRAM EQUDIFPC

parameter (MAXDATA=100)

real*8  tx(0:MAXDATA),ty(0:MAXDATA), x,xi,xf,yi
integer fi, n

    print *,' '
    write(*,"(' Input x begin: ')",advance='no')
    read *, xi
    write(*,"(' Input x end  : ')",advance='no')
    read *, xf
    write(*,"(' Input y begin: ')",advance='no')
    read *, yi
    write(*,"(' Input number of points: ')",advance='no')
    read *, n
    write(*,"(' Input finesse: ')",advance='no')
    read *, fi

    call equadiff_pc(tx,ty,xi,xf,yi,n,fi)

    !print results
    print *,' '
    print *,'      X                Y      '
    print *,' -----------------------------'
    do i=1, n 
      write(*,100)  tx(i),ty(i) 
    end do
    print *,' -----------------------------'
  
    stop
100 format(F12.6,'     ',F12.6)

end !of main program

  ! y'=(4xy+4x*sqrt(y))/(1+x²)
  real*8 Function f(x,y)
    real*8 x,y
    if (y>0.d0) then
      f = (4*x*y+4*x*dsqrt(y))/(1+x*x)
    else
      f = 0.d0
    end if
  end 

  !classical Runge-Kutta method of order 4
  Subroutine runge_kutta(h,x,y)
    real*8 h,x,y,a,b,c,d,f
    a=h*f(x,y)
    b=h*f(x+h/2.d0,y+a/2.d0)
    c=h*f(x+h/2.d0,y+b/2.d0)
    x=x+h
    d=h*f(x,y+c)
    y = y + (a+b+b+c+c+d)/6.d0
    return
  end

  !**********************************************************
  !*             Prediction-correction method               *
  !* ------------------------------------------------------ *
  !* INPUTS:                                                *
  !*	       xi      begin x value                        *
  !*	       xf      end x value                          *
  !*           y1      begin y value (at x=xi)              *
  !*           m       number of points to calculate        *
  !*           fi      finesse (number of intermediate      *
  !*                   points (for example 10)              *
  !* OUTPUTS:                                               *
  !*           tx      table of x values (m values)         *
  !*           ty      table of y values (m values)         * 
  !*                                                        *
  !* DESCRIPTION:                                           *
  !* The algorithm has the following steps:                 *
  !*     1. calculate y2, y3, y4 using a classical Runge-   *
  !*        Kutta method of order 4.                        *
  !*     2. from point (x4, y4), first estimate y(n+1) by   *
  !*        formula:                                        * 
  !*        y(n+1)=y(n) + H/24(55y'(n)-59y'(n-1)+37y'(n-2)  *
  !*               -9y'(n-3)                                *
  !*        then continue with formula:                     *
  !*        y(n+1)=y(n) + H/24(9y'(n+1)+19y'(n)-5y'(n-1)    *
  !*               +y'(n-2),                                *
  !*        noting that y'(n+1)=f(x(n+1),y(n+1)) with the   *
  !*        estimated value of y(n+1), until convergence is *
  !*        obtained.                                       *
  !**********************************************************
  Subroutine equadiff_pc(tx, ty, xi, xf, yi, m, fi)
    parameter(MAXDATA=100)
    real*8 tx(0:MAXDATA),ty(0:MAXDATA),xi,xf,yi,f
    integer m,fi
    integer*4 ni
    real*8 h,w,x,y,z,p(0:3)
    z=yi
    if (m>MAXDATA.or.fi<1) return
    h = (xf-xi)/fi/m
    p(3)=f(xi,yi)
    tx(0)=xi; ty(0)=yi
    k=0
    do i=1, m
      ni = (i-1)*fi-1
      do j=1, fi
	x=xi+h*dfloat(ni+j)
	k=k+1
	if (k<4) then
	  call runge_kutta(h,x,z)
	  p(3-k)=f(x,z)
	else 
	  x=x+h
          w=z+h/24*(55*p(0)-59*p(1)+37*p(2)-9*p(3))
10        continue
	  y=w
          w=z+h/24*(9*f(x,y)+19*p(0)-5*p(1)+p(2)) 
          if (dabs(y-w) > 1.d-10)  goto 10
	  z=w; p(3)=p(2); p(2)=p(1)
	  p(1)=p(0); p(0)=f(x,z)
        end if
      end do !j loop
      tx(i)=x; ty(i)=z
    end do !i loop
    return
  end


! end of file teqdifpc.f90