program main

c*********************************************************************72
c
cc fsolve_be_test() tests fsolve_be().
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

      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'fsolve_be_test():'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) '  Test fsolve_be()'
      write ( *, '(a)' ) '  which finds a vector ym that satisfies the'
      write ( *, '(a)' ) '  backward Euler residual:' 
      write ( *, '(a)' ) '    dydt(tm,ym) = (ym-yo)/(tm-to)'

      call predator_prey_test ( )
      call stiff_test ( )
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'fsolve_be_test():'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ' '
      call timestamp ( )

      stop 0
      end
      subroutine predator_prey_test ( )

c*********************************************************************72
c
cc predator_prey_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_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_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)

      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 stiff_test ( )

c*********************************************************************72
c
cc stiff_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_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_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 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