!**************************************************************************
!*  Solve a boundary value problem for a second order DE of the form:     *
!*  y" - p(x)*y = f(x) with conditions: y(x0)=a and y(x1)=b (x0<=x<=x1)   *
!*  using a Runge-Kutta integration method.                               *
!* ---------------------------------------------------------------------- *
!* Explanations: if y0(x) is a particular solution of y"-p(x)y = f(x) (1) *
!* and y1(x), y2(x) two lineary independant solutions of y"-p(x)y = 0 (2) *
!* (right hand side removed), then the general solution of (1) has the    *
!* form:  y(x) = y0(x) + c1*y1(x) + c2*y2(x), the constants c1, c2 are    *
!* obtained from boundary conditions.                                     *
!* A more efficient way to proceed is: 1) calculate a particular solution *
!* y1(x) of (2) such as y1(x0)=0 and 2) calculate a particular solution   *
!* y2(x) of (1) such as y2(x0)=a, then the general solution of (1) has the*
!* form:    y(x) = y1(x) + C * y2(x), constant C is determined by:        *
!*          C * y2(x1) + y1(x1) = b                                       *
!* ---------------------------------------------------------------------- *
!* SAMPLE RUN:                                                            *
!* Integrate y" - (1/x)y' - (3/x²)y = 5x² from x0=1 to x1=2 with boundary *
!* conditions: y(x0)=11 and y(x1)=36. The exact solution is:              *
!*           y(x) = 8/x + 2x² + x4                                        *
!* The initial 2nd order problem is replaced by the following first order *
!* DE systems:                                                            *
!*                 1st system                2nd system                   *
!*                 y' = z                    y' = z                       *
!*                 z' = z/x² + 3y/x²         z' = z/x² + 3y/x² + 5x²      *
!*                 with z'(x0)=0             with z'(x0)=11               *
!*                                                                        *
!*     X     Y estimated    Y true    Absolute Error                      *
!* --------------------------------------------------                     *
!*   1.0000  11.0000000   11.0000000     0.00E+00                         *
!*   1.0500  11.1498044   11.1498039     0.57E-06                         *
!*   1.1000  11.3988283   11.3988273     0.10E-05                         *
!*   1.1500  11.7472793   11.7472780     0.13E-05                         *
!*   1.2000  12.1962682   12.1962667     0.15E-05                         *
!*   1.2500  12.7476579   12.7476563     0.17E-05                         *
!*   1.3000  13.4039479   13.4039462     0.18E-05                         *
!*   1.3500  14.1681840   14.1681822     0.18E-05                         *
!*   1.4000  15.0438876   15.0438857     0.18E-05                         *
!*   1.4500  16.0349995   16.0349976     0.18E-05                         *
!*   1.5000  17.1458351   17.1458333     0.18E-05                         *
!*   1.5500  18.3810483   18.3810466     0.17E-05                         *
!*   1.6000  19.7456016   19.7456000     0.16E-05                         *
!*   1.6500  21.2447426   21.2447411     0.15E-05                         *
!*   1.7000  22.8839837   22.8839824     0.13E-05                         *
!*   1.7500  24.6690860   24.6690848     0.11E-05                         *
!*   1.8000  26.6060454   26.6060444     0.96E-06                         *
!*   1.8500  28.7010813   28.7010806     0.75E-06                         *
!*   1.9000  30.9606268   30.9606263     0.52E-06                         *
!*   1.9500  33.3913206   33.3913204     0.27E-06                         *
!*   2.0000  36.0000000   36.0000000     0.50E-13                         *
!*                                                                        *
!* ---------------------------------------------------------------------- *
!* REFERENCE: "Méthode de calcul numérique- Tome 2 - Programmes en Basic  *
!*             et en Pascal By Claude Nowakowski, Edition du P.S.I., 1984"*
!*             [BIBLI 04].                                                *
!*                                                                        *
!*                                    F90 Release By J-P Moreau, Paris.   *
!*                                           (www.jpmoreau.fr)            *
!**************************************************************************
Program Limits

integer, parameter :: SIZE = 100

real*8 X(0:SIZE), Y(0:3,0:SIZE),Z(0:3,0:SIZE)
real*8 cc,er,xl,x0,h,ya,yb,yex, FX
integer k,kl,l

  x0=1.d0         !starting x
  xl=2.d0         !ending x
  kl=20           !number of steps in x
  h=(xl-x0)/kl    !integration step
  ya=11.d0        !starting value for y
  yb=36.d0        !final value for y''
  X(0)=x0
  !integration loop
  do l=1, 2
    Y(l,0)=(l-1)*ya; Z(l,0)=1.d0
    do k=0, kl-1
      call RK4(l,X(k),Y(l,k),Z(l,k),h,X(k+1),Y(l,k+1),Z(l,k+1))
    end do
  end do
  !solution is y0=y2+cc*y1, cc is such as y(x1)=yb
  cc=(yb-Y(2,kl))/Y(1,kl)
  !write header
  print *,' ';
  print *,'     X     Y estimated    Y true    Absolute Error '
  print *,' --------------------------------------------------'
  !write kl+1 result lines
  do k=0, kl
    Y(0,k)=Y(2,k) + cc*Y(1,k)
    yex=FX(X(k))          !exact value
    er=DABS(yex-Y(0,k))   !current error
    Write(*,10) X(k), Y(0,k), yex, er
  end do
  print *,' --------------------------------------------------'
  stop
10 format(F9.4,F12.7,F14.7,E13.2)
  END

  !y'=z
  real*8 Function F(x,y,z)
  real*8 x,y,z
  F=z
  return
  End

  !l=1: z'=z/x²+3y/x² or l=2: z'=z/x²+3y/x²+5x²
  real*8 Function G(l,x,y,z)
  real*8 x,y,z,gg
  integer l
  gg=z/x+3.0*y/(x*x)
  if (l==2)  gg=gg+5.d0*x*x
  G=gg
  return
  End

  !exact solution: y=8/x + 2x3 + x4
  real*8 Function FX(x)
  real*8 x
  FX=8.d0/x+x*x*x*(x+2.d0)
  return
  End

  !Integrate sytem from x to x+h using Runge-Kutta
  Subroutine RK4(l,x,y,z,h,x1,y1,z1)
  real*8 x,y,z,h,x1,y1,z1
  real*8 c1,c2,c3,c4,d1,d2,d3,d4,h2
  real*8 F,G
  integer l
  c1=F(x,y,z)
  d1=G(l,x,y,z)
  h2=h/2.d0
  c2=F(x+h2,y+h2*c1,z+h2*d1)
  d2=G(l,x+h2,y+h2*c1,z+h2*d1)
  c3=F(x+h2,y+h2*c2,z+h2*d2)
  d3=G(l,x+h2,y+h2*c2,z+h2*d2)
  c4=F(x+h,y+h*c3,z+h*d3)
  d4=G(l,x+h,y+h*c3,z+h*d3)
  x1=x+h
  y1=y+h*(c1+2.d0*c2+2.d0*c3+c4)/6.d0
  z1=z+h*(d1+2.d0*d2+2.d0*d3+d4)/6.d0
  return
  End

!end of file limits.f90