subroutine stiff_deriv ( t, y, dydt )

!*****************************************************************************80
!
!! stiff_deriv() evaluates the right hand side of the stiff equation.
!
!  Discussion:
!
!    y' = lambda * ( cos(t-t0) - 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 ) t0
  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, t0_out = t0 )

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

  return
end
subroutine stiff_exact ( n, t, y )

!*****************************************************************************80
!
!! stiff_exact() evaluates the exact solution of the stiff ODE.
!
!  Discussion:
!
!    y' = lambda * ( cos(t-t0) - 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 ) mu
  real ( kind = rk8 ) t(n)
  real ( kind = rk8 ) t0
  real ( kind = rk8 ) y(n)
  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 )

  mu = y0 - lambda * lambda / ( lambda * lambda + 1.0 )

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

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