subroutine get_unit ( iunit )

!*****************************************************************************80
!
!! get_unit() returns a free FORTRAN unit number.
!
!  Discussion:
!
!    A "free" FORTRAN unit number is a value between 1 and 99 which
!    is not currently associated with an I/O device.  A free FORTRAN unit
!    number is needed in order to open a file with the OPEN command.
!
!    If IUNIT = 0, then no free FORTRAN unit could be found, although
!    all 99 units were checked (except for units 5, 6 and 9, which
!    are commonly reserved for console I/O).
!
!    Otherwise, IUNIT is a value between 1 and 99, representing a
!    free FORTRAN unit.  Note that GET_UNIT assumes that units 5 and 6
!    are special, and will never return those values.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    26 October 2008
!
!  Author:
!
!    John Burkardt
!
!  Output:
!
!    integer IUNIT, the free unit number.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer i
  integer ios
  integer iunit
  logical lopen

  iunit = 0

  do i = 1, 99

    if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then

      inquire ( unit = i, opened = lopen, iostat = ios )

      if ( ios == 0 ) then
        if ( .not. lopen ) then
          iunit = i
          return
        end if
      end if

    end if

  end do

  return
end
subroutine r8vec_linspace ( n, a, b, x )

!*****************************************************************************80
!
!! r8vec_linspace() creates a vector of linearly spaced values.
!
!  Discussion:
!
!    An R8VEC is a vector of R8's.
!
!    4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12.
!
!    In other words, the interval is divided into N-1 even subintervals,
!    and the endpoints of intervals are used as the points.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    14 March 2011
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of entries in the vector.
!
!    real ( kind = rk8 ) A, B, the first and last entries.
!
!  Output:
!
!    real ( kind = rk8 ) X(N), a vector of linearly spaced data.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) a
  real ( kind = rk8 ) b
  integer i
  real ( kind = rk8 ) x(n)

  if ( n == 1 ) then

    x(1) = ( a + b ) / 2.0D+00

  else

    do i = 1, n
      x(i) = ( real ( n - i,     kind = rk8 ) * a   &
             + real (     i - 1, kind = rk8 ) * b ) &
             / real ( n     - 1, kind = rk8 )
    end do

  end if

  return
end
subroutine stiff_deriv ( t, y, dydt )

!*****************************************************************************80
!
!! stiff_deriv() evaluates the right hand side of the stiff equation.
!
!  Discussion:
!
!    y' = lambda * ( cos(t) - y )
!    y(t0) = y0
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    23 June 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk8 ) T, Y: the time and solution value.
!
!  Output:
!
!    real ( kind = rk8 ) DYDT: the derivative value.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) dydt
  real ( kind = rk8 ) lambda
  real ( kind = rk8 ) t
  real ( kind = rk8 ) y

  interface
    subroutine stiff_parameters ( &
      lambda_in,  t0_in,  y0_in,  tstop_in, &
      lambda_out, t0_out, y0_out, tstop_out )
      integer, parameter :: rk8 = kind ( 1.0D+00 )
      real ( kind = rk8 ), optional :: lambda_in
      real ( kind = rk8 ), optional :: lambda_out
      real ( kind = rk8 ), optional :: t0_in
      real ( kind = rk8 ), optional :: t0_out
      real ( kind = rk8 ), optional :: y0_in
      real ( kind = rk8 ), optional :: y0_out
      real ( kind = rk8 ), optional :: tstop_in
      real ( kind = rk8 ), optional :: tstop_out
    end subroutine
  end interface

  call stiff_parameters ( lambda_out = lambda )

  dydt = lambda * ( cos ( t ) - y )

  return
end
subroutine stiff_euler ( n, t, y )

!*****************************************************************************80
!
!! stiff_euler() uses the Euler method on the stiff ODE.
!
!  Discussion:
!
!    y' = lambda * ( cos(t) - y )
!    y(t0) = y0
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    23 June 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N: the number of steps to take.
!
!  Output:
!
!    real ( kind = rk8 ) T(N+1), Y(N+1): the times and estimated solutions.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) dt
  integer i
  real ( kind = rk8 ) lambda
  real ( kind = rk8 ) t(n+1)
  real ( kind = rk8 ) t0
  real ( kind = rk8 ) tstop
  real ( kind = rk8 ) y(n+1)
  real ( kind = rk8 ) y0

  interface
    subroutine stiff_parameters ( &
      lambda_in,  t0_in,  y0_in,  tstop_in, &
      lambda_out, t0_out, y0_out, tstop_out )
      integer, parameter :: rk8 = kind ( 1.0D+00 )
      real ( kind = rk8 ), optional :: lambda_in
      real ( kind = rk8 ), optional :: lambda_out
      real ( kind = rk8 ), optional :: t0_in
      real ( kind = rk8 ), optional :: t0_out
      real ( kind = rk8 ), optional :: y0_in
      real ( kind = rk8 ), optional :: y0_out
      real ( kind = rk8 ), optional :: tstop_in
      real ( kind = rk8 ), optional :: tstop_out
    end subroutine
  end interface

  call stiff_parameters ( &
    lambda_out = lambda, t0_out = t0, y0_out = y0, tstop_out = tstop )

  dt = ( tstop - t0 ) / real ( n, kind = rk8 )

  t(1) = t0
  y(1) = y0

  do i = 1, n
    t(i+1) = t(i) + dt
    y(i+1) = y(i) + dt * lambda * ( cos ( t(i) ) - y(i) )
  end do

  return
end
subroutine stiff_euler_backward ( n, t, y )

!*****************************************************************************80
!
!! stiff_euler_backward() uses the backward Euler method on the stiff ODE.
!
!  Discussion:
!
!    y' = lambda * ( cos(t) - y )
!    y(t0) = y0
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    23 June 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N: the number of steps to take.
!
!  Output:
!
!    real ( kind = rk8 ) T(N+1), Y(N+1): the times and estimated solutions.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) dt
  integer i
  real ( kind = rk8 ) lambda
  real ( kind = rk8 ) t(n+1)
  real ( kind = rk8 ) t0
  real ( kind = rk8 ) tstop
  real ( kind = rk8 ) y(n+1)
  real ( kind = rk8 ) y0

  interface
    subroutine stiff_parameters ( &
      lambda_in,  t0_in,  y0_in,  tstop_in, &
      lambda_out, t0_out, y0_out, tstop_out )
      integer, parameter :: rk8 = kind ( 1.0D+00 )
      real ( kind = rk8 ), optional :: lambda_in
      real ( kind = rk8 ), optional :: lambda_out
      real ( kind = rk8 ), optional :: t0_in
      real ( kind = rk8 ), optional :: t0_out
      real ( kind = rk8 ), optional :: y0_in
      real ( kind = rk8 ), optional :: y0_out
      real ( kind = rk8 ), optional :: tstop_in
      real ( kind = rk8 ), optional :: tstop_out
    end subroutine
  end interface

  call stiff_parameters ( &
    lambda_out = lambda, t0_out = t0, y0_out = y0, tstop_out = tstop )

  dt = ( tstop - t0 ) / real ( n, kind = rk8 )

  t(1) = t0
  y(1) = y0

  do i = 1, n
    t(i+1) = t(i) + dt
    y(i+1) = ( y(i) + dt * lambda * cos ( t(i+1) ) ) &
      / ( 1.0D+00 + lambda * dt )
  end do

  return
end
subroutine stiff_exact ( n, t, y )

!*****************************************************************************80
!
!! stiff_exact() evaluates the exact solution of the stiff ODE.
!
!  Discussion:
!
!    y' = lambda * ( cos(t) - y )
!    y(t0) = y0
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    23 June 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N: the number of values.
!
!    real ( kind = rk8 ) T(N): the evaluation times.
!
!  Output:
!
!    real ( kind = rk8 ) Y(N): the exact solution values.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) lambda
  real ( kind = rk8 ) t(n)
  real ( kind = rk8 ) y(n)

  interface
    subroutine stiff_parameters ( &
      lambda_in,  t0_in,  y0_in,  tstop_in, &
      lambda_out, t0_out, y0_out, tstop_out )
      integer, parameter :: rk8 = kind ( 1.0D+00 )
      real ( kind = rk8 ), optional :: lambda_in
      real ( kind = rk8 ), optional :: lambda_out
      real ( kind = rk8 ), optional :: t0_in
      real ( kind = rk8 ), optional :: t0_out
      real ( kind = rk8 ), optional :: y0_in
      real ( kind = rk8 ), optional :: y0_out
      real ( kind = rk8 ), optional :: tstop_in
      real ( kind = rk8 ), optional :: tstop_out
    end subroutine
  end interface

  call stiff_parameters ( lambda_out = lambda )

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

  return
end
subroutine stiff_midpoint ( n, t, y )

!*****************************************************************************80
!
!! stiff_midpoint() uses the midpoint method on the stiff ODE.
!
!  Discussion:
!
!    y' = lambda * ( cos(t) - y )
!    y(t0) = y0
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    23 June 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N: the number of steps to take.
!
!  Output:
!
!    real ( kind = rk8 ) T(N+1), Y(N+1): the times and estimated solutions.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) dt
  integer i
  real ( kind = rk8 ) lambda
  real ( kind = rk8 ) t(n+1)
  real ( kind = rk8 ) t0
  real ( kind = rk8 ) tstop
  real ( kind = rk8 ) y(n+1)
  real ( kind = rk8 ) y0

  interface
    subroutine stiff_parameters ( &
      lambda_in,  t0_in,  y0_in,  tstop_in, &
      lambda_out, t0_out, y0_out, tstop_out )
      integer, parameter :: rk8 = kind ( 1.0D+00 )
      real ( kind = rk8 ), optional :: lambda_in
      real ( kind = rk8 ), optional :: lambda_out
      real ( kind = rk8 ), optional :: t0_in
      real ( kind = rk8 ), optional :: t0_out
      real ( kind = rk8 ), optional :: y0_in
      real ( kind = rk8 ), optional :: y0_out
      real ( kind = rk8 ), optional :: tstop_in
      real ( kind = rk8 ), optional :: tstop_out
    end subroutine
  end interface

  call stiff_parameters ( &
    lambda_out = lambda, t0_out = t0, y0_out = y0, tstop_out = tstop )

  dt = ( tstop - t0 ) / real ( n, kind = rk8 )

  t(1) = t0
  y(1) = y0

  do i = 1, n
    t(i+1) = t(i) + dt
    y(i+1) = &
      ( &
        y(i) + lambda * 0.5D+00 * dt * ( &
                                         cos ( t(i) ) - y(i) + cos ( t(i+1) ) &
                                       ) &
      ) &
      / ( 1.0D+00 + lambda * 0.5D+00 * dt )
  end do

  return
end
subroutine stiff_parameters ( &
  lambda_in,  t0_in,  y0_in,  tstop_in, &
  lambda_out, t0_out, y0_out, tstop_out )

!*****************************************************************************80
!
!! stiff_parameters() returns parameters of the stiff ODE.
!
!  Discussion:
!
!    y' = lambda * ( cos(t) - y )
!    y(t0) = y0
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    23 June 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk8 ) lambda_in: a parameter.
!
!    real ( kind = rk8 ) t0_in, y0_in: the initial time and value.
!
!    real ( kind = rk8 ) tstop_in: the final time.
!
!  Persistent:
!
!    real ( kind = rk8 ) lambda_default.
!
!    real ( kind = rk8 ) t0_default.
!
!    real ( kind = rk8 ) y0_default.
!
!    real ( kind = rk8 ) tstop_default.
!
!  Output:
!
!    real ( kind = rk8 ) lambda_out: a parameter.
!
!    real ( kind = rk8 ) t0_out, y0_out: the initial time and value.
!
!    real ( kind = rk8 ) tstop_out: the final time.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ), save     :: lambda_default = 50.0D+00
  real ( kind = rk8 ), optional :: lambda_in
  real ( kind = rk8 ), optional :: lambda_out
  real ( kind = rk8 ), save     :: t0_default = 0.0D+00
  real ( kind = rk8 ), optional :: t0_in
  real ( kind = rk8 ), optional :: t0_out
  real ( kind = rk8 ), save     :: y0_default = 0.0D+00
  real ( kind = rk8 ), optional :: y0_in
  real ( kind = rk8 ), optional :: y0_out
  real ( kind = rk8 ), save     :: tstop_default = 1.0D+00
  real ( kind = rk8 ), optional :: tstop_in
  real ( kind = rk8 ), optional :: tstop_out
!
!  New values, if supplied on input, overwrite the current values.
!
  if ( present ( lambda_in ) ) then
    lambda_default = lambda_in
  end if

  if ( present ( t0_in ) ) then
    t0_default = t0_in
  end if

  if ( present ( y0_in ) ) then
    y0_default = y0_in
  end if

  if ( present ( tstop_in ) ) then
    tstop_default = tstop_in
  end if
!
!  The current values are copied to the output.
!
  if ( present ( lambda_out ) ) then
    lambda_out = lambda_default
  end if

  if ( present ( t0_out ) ) then
    t0_out = t0_default
  end if

  if ( present ( y0_out ) ) then
    y0_out = y0_default
  end if

  if ( present ( tstop_out ) ) then
    tstop_out = tstop_default
  end if

  return
end