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 :: rk = 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 = rk ) A, B, the first and last entries.
!
!  Output:
!
!    real ( kind = rk ) X(N), a vector of linearly spaced data.
!
  implicit none

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

  integer n

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

  if ( n == 1 ) then

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

  else

    do i = 1, n
      x(i) = ( real ( n - i,     kind = rk ) * a   &
             + real (     i - 1, kind = rk ) * b ) &
             / real ( n     - 1, kind = rk )
    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:
!
!    06 April 2021
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk ) T, Y: the time and solution value.
!
!  Output:
!
!    real ( kind = rk ) DYDT: the derivative value.
!
  implicit none

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

  real ( kind = rk ) dydt
  real ( kind = rk ) lambda
  real ( kind = rk ) t
  real ( kind = rk ) t0
  real ( kind = rk ) y
  real ( kind = rk ) y0

  call stiff_parameters ( lambda, t0, y0 )

  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:
!
!    06 April 2021
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N: the number of steps to take.
!
!  Output:
!
!    real ( kind = rk ) T(N+1), Y(N+1): the times and estimated solutions.
!
  implicit none

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

  integer n

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

  call stiff_parameters ( lambda, t0, y0 )

  tstop = 1.0D+00

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

  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:
!
!    06 April 2021
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N: the number of steps to take.
!
!  Output:
!
!    real ( kind = rk ) T(N+1), Y(N+1): the times and estimated solutions.
!
  implicit none

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

  integer n

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

  call stiff_parameters ( lambda, t0, y0 )

  tstop = 1.0D+00

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

  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:
!
!    06 April 2021
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N: the number of values.
!
!    real ( kind = rk ) T(N): the evaluation times.
!
!  Output:
!
!    real ( kind = rk ) Y(N): the exact solution values.
!
  implicit none

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

  integer n

  real ( kind = rk ) lambda
  real ( kind = rk ) t(n)
  real ( kind = rk ) t0
  real ( kind = rk ) y(n)
  real ( kind = rk ) y0

  call stiff_parameters ( lambda, t0, y0 )

  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:
!
!    06 April 2021
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N: the number of steps to take.
!
!  Output:
!
!    real ( kind = rk ) T(N+1), Y(N+1): the times and estimated solutions.
!
  implicit none

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

  integer n

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

  call stiff_parameters ( lambda, t0, y0 )

  tstop = 1.0D+00

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

  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, t0, y0 )

!*****************************************************************************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:
!
!    06 April 2021
!
!  Author:
!
!    John Burkardt
!
!  Output:
!
!    real ( kind = rk ) LAMBDA, a parameter.
!
!    real ( kind = rk ) T0, Y0: the initial time and value.
!
  implicit none

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

  real ( kind = rk ) lambda
  real ( kind = rk ) t0
  real ( kind = rk ) y0

  lambda = 50.0D+00
  t0 = 0.0D+00
  y0 = 0.0D+00

  return
end