!********************************************************************
!*             Differential equations of order N                    *
!*             by Runge-Kutta method of order 4                     *
!*                with boundary conditions                          *
!* ---------------------------------------------------------------- *
!* Reference: 'Analyse en Turbo Pascal versions 5.5 et 6.0' by Marc *
!*             DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991.    *
!*                                                                  *
!*                             F90 Version By J-P Moreau, Paris     *
!*                                     (www.jpmoreau.fr)            *
!* ---------------------------------------------------------------- *
!* SAMPLE RUN:                                                      *
!*                                                                  *
!* Example: integrate y"=(2y'y'+yy)/y from x=4 to x=6               *
!*          with initial conditions: y(4)=2 and y'(4)=-2            *
!*          and final condition y(6)=2                              *
!*                                                                  *
!* Exact solution is:   y = 2cos(1)/cos(x-5)                        *
!*                                                                  *
!*        DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER N          *
!*              of type y(n) = f(y(n-1),...,y',y,x)                 *
!*                 with boundary conditions                         *
!*                                                                  *
!*   Order:             2                                           *
!*   Begin value x    : 4.00000000000000                            *
!*   end value x      : 6.00000000000000                            *
!*   Y value at begin  x   :  2.00000000000000                      *
!*   Y' value at begin x   : -2.00000000000000                      *
!*   Desired value at end x:  2.00000000000000                      *   
!*   number of points   :      11                                   *
!*   finesse            :      20                                   *
!*   Max. iterations    :     100                                   *
!*                                                                  *
!*   Last iteration:                                                *
!*   u  =-3.11481552124023        du  =-7.629394531250000E-007      *
!*   dyb= 8.195419964884820E-008           67 steps                 *
!*                                                                  *
!*       X            Y                                             *
!* ---------------------------                                      *
!*   4.000000     2.000000                                          *
!*   4.200000     1.551018                                          *
!*   4.400000     1.309291                                          *
!*   4.600000     1.173217                                          *
!*   4.800000     1.102583                                          *
!*   5.000000     1.080605                                          *
!*   5.200000     1.102583                                          *
!*   5.400000     1.173217                                          *
!*   5.600000     1.309291                                          *
!*   5.800000     1.551018                                          *
!*   6.000000     2.000000                                          *
!*                                                                  *
!********************************************************************
PROGRAM TSHOOT

integer fi
real*8  vect(1:50), xi,xf,yi(10)

real*8  du,dyb,olddyb,u,yb

  print *,' '
  print *,'      DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER N'
  print *,'            of type y(n) = f(y(n-1),...,y'',y,x)'
  print *,'               with boundary conditions'
  print *,' '
  iordre=2
  print *,' Order: ', iordre
  xi=4.d0; xf=6.d0
  print *,' Begin value x: ', xi
  print *,' End value x  : ', xf
  yi(0)=2.d0; yi(1)=-2.d0
  print *,' Y value at begin  x   :', yi(0)
  print *,' Y'' value at begin x   :', yi(1)
  yb=2.d0
  print *,' Desired value at end x:', yb
  ndata=11
  print *,' Number of points:', ndata 
  fi=20
  print *,' Finesse         :', fi
  itermax=100
  print *,' Max. iterations :', itermax

  vect=0.d0

  u=yi(1)          !First derivative initial value
  du=dabs(u)/10.d0
  niter=1
! estimate initial error on desired final value
  call Equadifn(vect,xi,xf,yi,ndata,iordre,fi,.false.)

  olddyb=dabs(vect(ndata)-yb)
  dyb=olddyb

! Iterate until desired final value yb is obtained
  do while(dyb>1.d-7.and.niter<=itermax)
    u=u+du; yi(1)=u
    call Equadifn(vect,xi,xf,yi,ndata,iordre,fi,.false.)
    dyb=dabs(vect(ndata)-yb)
    if (dabs(dyb)>dabs(olddyb)) then  !overshoot?
      if (niter==1) then
        du=-du
      else
        du=-du/2.d0
      end if
    end if
    olddyb=dyb 
	niter=niter+1
  end do

  print *,' '
  print *,' Last iteration:'
  print *,' u  =',u,'  du  =',du
  print *,' dyb=',dyb,'  ',niter,' steps'

  !call subroutine equadifn() and print results
  call Equadifn(vect,xi,xf,yi,ndata,iordre,fi,.true.)

  stop

50 format('    y',i1,' value at x0   : ')

END

!Example: y'=(2y'y'+yy)/y
real*8 FUNCTION fp(x,y)
real*8 x,y(10)
  if (dabs(y(0))<1.d-12)  y(0)=1.d-12
  fp=(2.d0*y(1)*y(1)+y(0)*y(0))/y(0)
end

!***************************************************************************
!*        SOLVING DIFFERENTIAL EQUATIONS WITH 1 VARIABLE OF ORDER N        *
!*                of type y(n) = f(y(n-1),y(n-2),...,y',y,x)               *
!* ----------------------------------------------------------------------- *
!*  INPUTS:                                                                *
!*    xi, xf    Begin, end values of variable x                            *
!*    yi        Begin values at xi (of f and derivatives)                  *
!*    m         Number of points to calculate                              *
!*    n         Order of differential equation                             *
!*    fi        finesse (number of intermediary points)                    *
!*                                                                         *
!*  OUTPUTS:                                                               *
!*    t         real vector storing m results for function y               *
!* ----------------------------------------------------------------------- *
!*  EXAMPLE:    y' = (2 y'y' + yy) / y with y(4) = 2, y'(4) = -2*tan(1)    *
!*              Exact solution:  y = (2 cos(1))/cos(x-5)                   *
!***************************************************************************
Subroutine Equadifn(t,xi,xf,yi,m,n,fi,flag)
integer fi
logical flag
real*8  t(*),xi,xf,yi(*), fp
real*8  h,x,a,b,c,d
real*8  ta(10),tb(10),tc(10),td(10),y(10),z(10)

!  print *,' ** m=', m
!  print *,' ** ordre=', n
!  print *,' ** fi=', fi

   if (fi<1) return
   h = (xf - xi) / fi / dfloat(m-1)
   n1=n-1; n2=n-2
   t(1)=yi(0)
   do k=0, n1
     y(k)=yi(k); z(k)=yi(k)
   end do
   do i=1, m
     ni=(i-1)*fi-1
     do j=1, fi
       x=xi+h*(ni+j)
       do k=0, n1
         y(k)=z(k)
       end do
       a=h*fp(x,y)
       do k=0, n2
         ta(k)=h*y(k+1); y(k)=z(k)+ta(k)/2.d0
       end do
       y(n1)=z(n1)+a/2.d0
       x=x+h/2
       b=h*fp(x,y)
       do k=0, n2
         tb(k)=h*y(k+1); y(k)=z(k)+tb(k)/2.d0
       end do
       y(n1)=z(n1)+b/2.d0
       c=h*fp(x,y)
       do k=0, n2
         tc(k)=h*y(k+1); y(k)=z(k)+tc(k)
       end do
       y(n1)=z(n1)+c
       x=x+h/2
       d=h*fp(x,y)
       do k=0, n2
         z(k)=z(k)+(ta(k)+2.d0*tb(k)+2.d0*tc(k)+h*y(k+1))/6.d0
       end do
       z(n1)=z(n1)+(a+b+b+c+c+d)/6.d0
     end do
     t(i+1)=z(0)
   end do
   !call subroutine Display if required
   if (flag) call Display(m,xi,xf,t)
   return
end

!print solutions
Subroutine Display(m,xi,xf,t)
  integer i,m
  real*8 h,x,xi,xf,t(*)
  h=(xf-xi)/dfloat(m-1)
  x=xi-h
  print *,' '
  print *, '      X         Y     '
  print *, '----------------------'
  do i=1, m
    x=x+h
    write(*,100)  x, t(i)
  end do
  print *,' '
  return
100 format(' ',f9.6,'  ',f9.6)  
end

!End of file tshoot.f90