!>
!> @file difequ.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
subroutine difequ ( mode, xx, v )

!*************************************************************************
!
!! DIFEQU returns information about a differential equation.
!
!  Reference:
!
!    Carl DeBoor,
!    A Practical Guide to Splines,
!    Springer Verlag.
!
!  Parameters:
!
!    Input, integer MODE, an integer indicating the task to be performed.
!    1, initialization
!    2, evaluate  de  at  xx
!    3, specify the next side condition
!    4, analyze the approximation
!
!    Input, real ( kind = 8 ) XX, a point at which information is wanted
!
!    Output, real ( kind = 8 ) V, depends on the MODE.  See comments below
!
  implicit none

  integer, parameter :: npiece = 100
  integer, parameter :: ncoef = 2000

  real ( kind = 8 ) break
  real ( kind = 8 ) coef
  real ( kind = 8 ) eps
  real ( kind = 8 ) ep1
  real ( kind = 8 ) ep2
  real ( kind = 8 ) error
  real ( kind = 8 ) factor
  integer i
  integer iside
  integer itermx
  integer k
  integer kpm
  integer l
  integer m
  integer mode
  real ( kind = 8 ) rho
  real ( kind = 8 ) s2ovep
  real ( kind = 8 ) solutn
  real ( kind = 8 ) un
  real ( kind = 8 ) v(20)
  real ( kind = 8 ) value
  real ( kind = 8 ) x
  real ( kind = 8 ) xside
  real ( kind = 8 ) xx

  common /approx/ break(npiece),coef(ncoef),l,kpm
  common /side/ m,iside,xside(10)
  common /other/ itermx,k,rho(19)
!
!  This sample of DIFEQU is for the example in chapter xv.  It is a
!  nonlinear second order two point boundary value problem.
!
  go to (10,50,60,110), mode
!
!  Initialize everything,  Set the order M of the differential equation,
!  the nondecreasing sequence xside(i),i=1,...,m, of points at which side
!  conditions are given and anything else necessary.
!
   10 continue

  m = 2
  xside(1) = 0.0D+00
  xside(2) = 1.0D+00
!
!  Print out heading.
!
  write ( *, * ) ' '
  write ( *, * ) '  Carrier''s nonlinear perturb. problem'
  write ( *, * ) ' '

  eps = 0.005D+00
  write(*,*)'EPS = ',eps
!
!  Set constants used in formula for solution below.
!
  factor = ( sqrt ( 2.0D+00 ) + sqrt ( 3.0D+00 ) )**2
  s2ovep = sqrt ( 2.0D+00 / eps )
!
!  Initial guess for Newton iteration. un(x)=x*x-1.
!
  l = 1
  break(1) = 0.0D+00
  do i = 1, kpm
    coef(i) = 0.0D+00
  end do
  coef(1) = -1.0D+00
  coef(3) = 2.0D+00
  itermx = 10
  return
!
!  Provide value of left side coefficients and right side at  xx .
!  specifically, at  xx  the dif.equ. reads:
!
!    v(m+1)d**m+v(m)d**(m-1) + ... + v(1)d**0 = v(m+2)
!
!  in terms of the quantities v(i),i=1,...,m+2, to be computed here.
!
   50 continue

  v(3) = eps
  v(2) = 0.0D+00
  call ppvalu(break,coef,l,kpm,xx,0,un)
  v(1) = 2.0D+00 * un
  v(4) = un**2 + 1.0D+00
  return
!
!  provide the M side conditions. these conditions are of the form
!        v(m+1)d**m+v(m)d**(m-1) + ... + v(1)d**0 = v(m+2)
!  in terms of the quantities v(i),i=1,...,m+2, to be specified here.
!  note that v(m+1)=0  for customary side conditions.
!
   60 continue

  v(m+1) = 0.0D+00
  if ( iside == 1 ) then
    v(2) = 1.0D+00
    v(1) = 0.0D+00
    v(4) = 0.0D+00
    iside = iside+1
  else if ( iside == 2 ) then
    v(2) = 0.0D+00
    v(1) = 1.0D+00
    v(4) = 0.0D+00
    iside = iside + 1
  end if

  return
!
!  Calculate the error near the boundary layer at 1.
!
  110 continue

  write(*,*)' '
  write(*,*)' X, G(X) and G(X)-F(X) at selected points:'
  write(*,*)' '

  x = 0.75D+00

  do i = 1, 9
    ep1 = exp ( s2ovep * ( 1.0D+00 - x ) ) * factor
    ep2 = exp ( s2ovep * ( 1.0D+00 + x ) ) * factor
    solutn = 12.0D+00 / ( 1.0D+00 + ep1 )**2 * ep1 &
            +12.0D+00 / ( 1.0D+00 + ep2 )**2 * ep2 - 1.0D+00
    call ppvalu(break,coef,l,kpm,x,0,value)
    error = solutn-value
    write ( *, '(1x,3g14.6)' ) x, solutn, error
    x = x+0.03125
  end do

  return
end