program main

c*********************************************************************72
c
cc fsolve_test() tests fsolve().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    30 November 2023
c
c  Author:
c
c    John Burkardt
c
      implicit none

      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'fsolve_test():'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) '  Test fsolve()'
      write ( *, '(a)' ) '  solve systems of nonlinear equations.'

      call fsolve_test1 ( )
      call fsolve_test2 ( )
      call fsolve_test3 ( )
      call fsolve_test4 ( )

      call predator_prey_be_test ( )
      call predator_prey_tr_test ( )

      call stiff_bdf2_test ( )
      call stiff_be_test ( )
      call stiff_tr_test ( )
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'fsolve_test():'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ' '
      call timestamp ( )

      stop 0
      end
      subroutine fsolve_test1 ( )

c*********************************************************************72
c
cc fsolve_test1() tests fsolve() on a system of 1 equation.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 October 2023
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer, parameter :: n = 1

      external f1
      double precision fx(n)
      integer info
      double precision tol
      double precision x(n)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'fsolve_test1():'
      write ( *, '(a)' ) '  Solves a nonlinear system of 1 equation.'

      x(1:1) = (/ 0.0D+00 /)

      call f1 ( n, x, fx )      
      call r8vec2_print ( n, x, fx, '  Initial X and F(X)' )

      tol = 0.00001D+00

      call fsolve ( f1, n, x, fx, tol, info )

      write ( *, '(a)' ) ' '
      if ( info == 1 ) then
        write ( *, '(a)' ) '  Satisfactory computation.'
      else
        write ( *, '(a,i6)' ) '  Returned value of INFO = ', info
      end if

      call r8vec2_print ( n, x, fx, '  Final X and F(X)' )

      return
      end
      subroutine f1 ( n, x, fx )

c*********************************************************************72
c
cc f1() evaluates a nonlinear system of 1 equation.
c
c  Discussion:
c
c    This is Kepler's equation.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 October 2023
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer N, the number of variables.
c
c    double precision X(N), the variable values.
c
c  Output:
c
c    double precision FX(N), the function values at X.
c
      implicit none

      integer n

      double precision e
      double precision fx(n)
      double precision m
      double precision x(n)

      e = 0.8D+00
      m = 5.0D+00
      fx(1) = x(1) - m - e * sin ( x(1) )

      return
      end
      subroutine fsolve_test2 ( )

c*********************************************************************72
c
cc fsolve_test2() tests fsolve() on a system of 2 equations.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 October 2023
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer, parameter :: n = 2

      external f2
      double precision fx(n)
      integer info
      double precision tol
      double precision x(n)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'fsolve_test2():'
      write ( *, '(a)' ) '  Solves a nonlinear system of 2 equations.'

      x(1:2) = (/ 3.0D+00, 0.0D+00 /)

      call f2 ( n, x, fx )      
      call r8vec2_print ( n, x, fx, '  Initial X and F(X)' )

      tol = 0.00001D+00

      call fsolve ( f2, n, x, fx, tol, info )

      write ( *, '(a)' ) ' '
      if ( info == 1 ) then
        write ( *, '(a)' ) '  Satisfactory computation.'
      else
        write ( *, '(a,i6)' ) '  Returned value of INFO = ', info
      end if

      call r8vec2_print ( n, x, fx, '  Final X and F(X)' )

      return
      end
      subroutine f2 ( n, x, fx )

c*********************************************************************72
c
cc f2() evaluates a nonlinear system of 2 equations.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 October 2023
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer N, the number of variables.
c
c    double precision X(N), the variable values.
c
c  Output:
c
c    double precision FX(N), the function values at X.
c
      implicit none

      integer n

      double precision fx(n)
      double precision x(n)

      fx(1) = x(1) * x(1) - 10.0D+00 * x(1) + x(2) * x(2) + 8.0D+00
      fx(2) = x(1) * x(2) * x(2) + x(1) - 10.0D+00 * x(2) + 8.0D+00

      return
      end
      subroutine fsolve_test3 ( )

c*********************************************************************72
c
cc fsolve_test3() tests fsolve() on a system of 4 equations.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 October 2023
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer, parameter :: n = 4

      external f3
      double precision fx(n)
      integer info
      double precision tol
      double precision x(n)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'fsolve_test3():'
      write ( *, '(a)' ) '  Solves a nonlinear system of 4 equations.'

      x(1:n) = 0.0D+00

      call f3 ( n, x, fx )      
      call r8vec2_print ( n, x, fx, '  Initial X and F(X)' )

      tol = 0.00001D+00

      call fsolve ( f3, n, x, fx, tol, info )

      write ( *, '(a)' ) ' '
      if ( info == 1 ) then
        write ( *, '(a)' ) '  Satisfactory computation.'
      else
        write ( *, '(a,i6)' ) '  Returned value of INFO = ', info
      end if

      call r8vec2_print ( n, x, fx, '  Final X and F(X)' )

      return
      end
      subroutine f3 ( n, x, fx )

c*********************************************************************72
c
cc f3() evaluates a nonlinear system of 4 equations.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 October 2023
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer N, the number of variables.
c
c    double precision X(N), the variable values.
c
c  Output:
c
c    double precision FX(N), the function values at X.
c
      implicit none

      integer n

      double precision fx(n)
      integer i
      double precision x(n)

      do i = 1, n
        fx(i) = ( x(i) - i )**2
      end do

      return
      end
      subroutine fsolve_test4 ( )

c*********************************************************************72
c
cc fsolve_test4() tests fsolve() on a system of 8 equations.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 October 2023
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer, parameter :: n = 8

      external f4
      double precision fx(n)
      integer info
      double precision tol
      double precision x(n)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'fsolve_test4():'
      write ( *, '(a)' ) '  Solve a nonlinear system of 8 equations.'

      x(1:n) = 0.0D+00

      call f4 ( n, x, fx )      
      call r8vec2_print ( n, x, fx, '  Initial X and F(X)' )

      tol = 0.00001D+00

      call fsolve ( f4, n, x, fx, tol, info )

      write ( *, '(a)' ) ' '
      if ( info == 1 ) then
        write ( *, '(a)' ) '  Satisfactory computation.'
      else
        write ( *, '(a,i6)' ) '  Returned value of INFO = ', info
      end if

      call r8vec2_print ( n, x, fx, '  Final X and F(X)' )

      return
      end
      subroutine f4 ( n, x, fx )

c*********************************************************************72
c
cc f4() evaluates a nonlinear system of 8 equations.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 October 2023
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer N, the number of variables.
c
c    double precision X(N), the variable values.
c
c  Output:
c
c    double precision FX(N), the function values at X.
c
      implicit none

      integer n

      double precision fx(n)
      integer i
      double precision x(n)

      do i = 1, n

        fx(i) = ( 3.0D+00 - 2.0D+00 * x(i) ) * x(i) + 1.0D+00

        if ( 1 < i ) then
          fx(i) = fx(i) - x(i-1)
        end if

        if ( i < n ) then
          fx(i) = fx(i) - 2.0D+00 * x(i+1)
        end if

      end do

      return
      end
      subroutine predator_prey_be_test ( )

c*********************************************************************72
c
cc predator_prey_be_test() tests fsolve_be() on the predator prey ODE.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    10 November 2023
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer, parameter :: m = 2

      external predator_prey_dydt
      double precision fm(m)
      integer info
      double precision tm
      double precision to
      double precision tol
      double precision ym(m)
      double precision yo(m)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'predator_prey_be_test():'

      to = 0.0D+00
      yo = (/ 5000.0D+00, 100.0D+00 /)
      tm = 1.0D+00
      ym = yo
      tol = 1.0D-05

      call backward_euler_residual ( predator_prey_dydt, m, to, yo, 
     &  tm, ym, fm )
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Initial ||ym-yo-(tm-to)*dydt(tm,ym)||:'
      write ( *, '(g14.6)' ) sqrt ( sum ( fm**2 ) )

      call fsolve_be ( predator_prey_dydt, m, to, yo, tm, ym, fm, 
     &  tol, info )
      if ( info /= 1 ) then
        write ( *, '(a)' ) ''
        write ( *, '(a,i6)' ) 
     &    '  fsolve_be() returned error flag info = ', info
      end if

      call backward_euler_residual ( predator_prey_dydt, m, to, yo, 
     &  tm, ym, fm )
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Final ||ym-yo-(tm-to)*dydt(tm,ym)||:'
      write ( *, '(g14.6)' ) sqrt ( sum ( fm**2 ) )

      return
      end
      subroutine predator_prey_tr_test ( )

c*********************************************************************72
c
cc predator_prey_tr_test() tests fsolve_tr() on the predator prey ODE.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    30 November 2023
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer, parameter :: m = 2

      external predator_prey_dydt
      double precision fm(m)
      integer info
      double precision tn
      double precision to
      double precision tol
      double precision yn(m)
      double precision yo(m)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'predator_prey_tr_test():'

      to = 0.0D+00
      yo = (/ 5000.0D+00, 100.0D+00 /)
      tn = 1.0D+00
      yn = yo
      tol = 1.0D-05

      call trapezoidal_residual ( predator_prey_dydt, m, to, yo, 
     &  tn, yn, fm )
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Initial residual:'
      write ( *, '(g14.6)' ) sqrt ( sum ( fm**2 ) )

      call fsolve_tr ( predator_prey_dydt, m, to, yo, tn, yn, fm, 
     &  tol, info )
      if ( info /= 1 ) then
        write ( *, '(a)' ) ''
        write ( *, '(a,i6)' ) 
     &    '  fsolve_tr() returned error flag info = ', info
      end if

      call trapezoidal_residual ( predator_prey_dydt, m, to, yo, 
     &  tn, yn, fm )
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Final residual:'
      write ( *, '(g14.6)' ) sqrt ( sum ( fm**2 ) )

      return
      end
      subroutine predator_prey_dydt ( t, y, dydt )

c*********************************************************************72
c
cc predator_prey_dydt() evaluates the right hand side of the system.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    22 February 2020
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    George Lindfield, John Penny,
c    Numerical Methods Using MATLAB,
c    Second Edition,
c    Prentice Hall, 1999,
c    ISBN: 0-13-012641-1,
c    LC: QA297.P45.
c
c  Input:
c
c    double precision t: the current time.
c
c    double precision y(2): the current solution variables, rabbits and foxes.
c
c  Output:
c
c    double precision dydt(2): the right hand side of the 2 ODE's.
c
      implicit none

      double precision dydt(2)
      double precision t
      double precision y(2)

      call r8_fake_use ( t )

      dydt(1) =    2.0 * y(1) - 0.001 * y(1) * y(2)
      dydt(2) = - 10.0 * y(2) + 0.002 * y(1) * y(2)

      return
      end
      subroutine r8_fake_use ( x )

c*****************************************************************************80
c
cc r8_fake_use() pretends to use an R8 variable.
c
c  Discussion:
c
c    Some compilers will issue a warning if a variable is unused.
c    Sometimes there's a good reason to include a variable in a program,
c    but not to use it.  Calling this function with that variable as
c    the argument will shut the compiler up.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    19 August 2023
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    double precision x: the variable to be "used".
c
      implicit none

      double precision x

      if ( x /= x ) then
        write ( *, '(a)' ) '  r8_fake_use(): variable is NAN.'
      end if

      return
      end
      subroutine r8vec2_print ( n, a1, a2, title )

c*********************************************************************72
c
cc r8vec2_print() prints an R8VEC2.
c
c  Discussion:
c
c    An R8VEC2 is a dataset consisting of N pairs of R8s, stored
c    as two separate vectors A1 and A2.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    06 February 2008
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, integer N, the number of components of the vector.
c
c    Input, double precision A1(N), A2(N), the vectors to be printed.
c
c    Input, character * ( * ) TITLE, a title.
c
      implicit none

      integer n

      double precision a1(n)
      double precision a2(n)
      integer i
      character * ( * ) title

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) trim ( title )
      write ( *, '(a)' ) ' '
      do i = 1, n
        write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6)' ) i, ':', a1(i), a2(i)
      end do

      return
      end
      subroutine stiff_bdf2_test ( )

c*********************************************************************72
c
cc stiff_bdf2_test() tests fsolve_bdf2() on the stiff ODE.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    30 November 2023
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer, parameter :: m = 1

      external stiff_dydt
      double precision fm(m)
      integer info
      double precision t1
      double precision t2
      double precision t3
      double precision tol
      double precision y1(m)
      double precision y2(m)
      double precision y3(m)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'stiff_bdf2_test():'

      t1 = 0.0D+00
      y1(1) = 0.0D+00
      t2 = 0.25
      y2(1) = 50.0D+00 * ( sin ( t2 ) + 50.0D+00 * cos ( t2 ) 
     &  - 50.0D+00 * exp ( - 50.0D+00 * t2 ) ) / 2501.0D+00
      t3 = 0.5
      y3 = y2
      tol = 1.0D-05

      call bdf2_residual ( stiff_dydt, m, t1, y1, 
     &  t2, y2, t3, y3, fm )
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Initial residual:'
      write ( *, '(g14.6)' ) sqrt ( sum ( fm**2 ) )

      call fsolve_bdf2 ( stiff_dydt, m, t1, y1, t2, y2, t3, y3,
     &  fm, tol, info )

      if ( info /= 1 ) then
        write ( *, '(a)' ) ''
        write ( *, '(a,i6)' ) 
     &    '  fsolve_bdf2() returned error flag info = ', info
      end if

      call bdf2_residual ( stiff_dydt, m, t1, y1, 
     &  t2, y2, t3, y3, fm )
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Final residual:'
      write ( *, '(g14.6)' ) sqrt ( sum ( fm**2 ) )

      return
      end
      subroutine stiff_be_test ( )

c*********************************************************************72
c
cc stiff_be_test() tests fsolve_be() on the stiff ODE.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    08 November 2023
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer, parameter :: m = 1

      external stiff_dydt
      double precision fm(m)
      integer info
      double precision tm
      double precision to
      double precision tol
      double precision ym(m)
      double precision yo(m)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'stiff_be_test():'

      to = 0.0D+00
      yo = (/ 0.0D+00 /)
      tm = 0.25
      ym = yo
      tol = 1.0D-05

      call backward_euler_residual ( stiff_dydt, m, to, yo, 
     &  tm, ym, fm )
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Initial ||ym-yo-(tm-to)*dydt(tm,ym)||:'
      write ( *, '(g14.6)' ) sqrt ( sum ( fm**2 ) )

      call fsolve_be ( stiff_dydt, m, to, yo, tm, ym, fm, tol, info )

      if ( info /= 1 ) then
        write ( *, '(a)' ) ''
        write ( *, '(a,i6)' ) 
     &    '  fsolve_be() returned error flag info = ', info
      end if

      call backward_euler_residual ( stiff_dydt, m, to, yo, 
     &  tm, ym, fm )
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Final ||ym-yo-(tm-to)*dydt(tm,ym)||:'
      write ( *, '(g14.6)' ) sqrt ( sum ( fm**2 ) )

      return
      end
      subroutine stiff_tr_test ( )

c*********************************************************************72
c
cc stiff_tr_test() tests fsolve_tr() on the stiff ODE.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    30 November 2023
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer, parameter :: m = 1

      external stiff_dydt
      double precision fm(m)
      integer info
      double precision tn
      double precision to
      double precision tol
      double precision yn(m)
      double precision yo(m)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'stiff_tr_test():'

      to = 0.0D+00
      yo(1) = 0.0D+00
      tn = 0.25
      yn(1) = yo(1)
      tol = 1.0D-05

      call trapezoidal_residual ( stiff_dydt, m, to, yo, 
     &  tn, yn, fm )
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Initial residual:'
      write ( *, '(g14.6)' ) sqrt ( sum ( fm**2 ) )

      call fsolve_tr ( stiff_dydt, m, to, yo, tn, yn,
     &  fm, tol, info )

      if ( info /= 1 ) then
        write ( *, '(a)' ) ''
        write ( *, '(a,i6)' ) 
     &    '  fsolve_tr() returned error flag info = ', info
      end if

      call trapezoidal_residual ( stiff_dydt, m, to, yo, 
     &  tn, yn, fm )
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Final residual:'
      write ( *, '(g14.6)' ) sqrt ( sum ( fm**2 ) )

      return
      end
      subroutine stiff_dydt ( t, y, dydt )

c*********************************************************************72
c
cc stiff_dydt() evaluates the right hand side of the stiff ODE.
c
c  Discussion:
c
c    y' = 50 * ( cos(t) - y )
c    y(0) = 0
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    10 November 2023
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    double precision T, Y: the time and solution value.
c
c  Output:
c
c    double precision DYDT: the derivative value.
c
      implicit none

      double precision dydt(1)
      double precision t
      double precision y(1)

      dydt(1) = 50.0D+00 * ( cos ( t ) - y(1) )

      return
      end
      subroutine stiff_exact ( n, t, y )

c*****************************************************************************80
c
cc stiff_exact() evaluates the exact solution of the stiff ODE.
c
c  Discussion:
c
c    y' = 50 * ( cos(t) - y )
c    y(0) = 0
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    27 April 2020
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer N: the number of values.
c
c    double precision T(N): the evaluation times.
c
c  Output:
c
c     double precision Y(N): the exact solution values.
c
      implicit none

      integer n

      double precision t(n)
      double precision y(n)

      y(1:n) = 50.0D+00 * ( sin ( t ) + 50.0D+00 * cos(t) 
     &  - 50.0D+00 * exp ( - 50.0D+00 * t ) ) / 2501.0D+00

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc timestamp() prints the YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 June 2014
c
c  Author:
c
c    John Burkardt
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *, 
     &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) 
     &  d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, 
     &  trim ( ampm )

      return
      end