#! /usr/bin/env python3
#
def stiff_ode ( x, y ):

#*****************************************************************************80
#
## stiff_ode() evaluates the right hand side of the stiff ODE.
#
#  input:
#    x is independent variable
#    y is dependent variable
#  output:
#    dydx is the value of f_ode(x,y).
#
  import numpy as np
  LAMBDA = stiff_lambda ( )
  dydx = LAMBDA * ( - y + np.sin ( x ) ) 
  return dydx

def stiff_solution ( x ):

#*****************************************************************************80
#
## stiff_solution() evaluates the exact solution of the stiff ODE.
#
#  input:
#    x is the independent variable
#  output:
#    y is the solution value
#
  import numpy as np
  LAMBDA = stiff_lambda ( )
  y = ( LAMBDA**2 / ( 1 + LAMBDA**2 ) ) * np.sin(x) + \
      ( LAMBDA    / ( 1 + LAMBDA**2 ) ) * ( np.exp ( - LAMBDA * x ) - np.cos ( x ) )
  return y

def stiff_lambda ( lambda_input = None ):

#*****************************************************************************80
#
## stiff_lambda() allows the user to work with a global variable lambda.
#
#  Discussion:
#
#    Don't actually call the variable "lambda" because that's a Python keyword!
#
#  Modified:
#
#    02 October 2021
#
#  Author:
#
#    John Burkardt
#
  if ( not hasattr ( stiff_lambda, "lambda_default" ) ):
    stiff_lambda.lambda_default = 4.0

  if ( lambda_input is not None ):
    stiff_lambda.lambda_default = lambda_input

  lambda_output = stiff_lambda.lambda_default

  return lambda_output

def stiff_lambda_test ( ):

#*****************************************************************************80
#
## stiff_lambda_test() demonstrates the use of stiff_lambda().
#
  print ( "stiff_lambda_test:" )
  
  lambda_output = stiff_lambda ( )
  print ( "" )
  print ( "  Calling with no arguments returns the default:" )
  print ( "    lambda_output = stiff_lambda ( ) = ", lambda_output )

  stiff_lambda ( 99.0 )
  print ( "" )
  print ( "  Calling with an argument resets the default:" )
  print ( "    stiff_lambda ( 99.0 )" )

  lambda_output = stiff_lambda ( )
  print ( "" )
  print ( "  Calling with no arguments returns the (new) default:" )
  print ( "    lambda_output = stiff_lambda ( ) = ", lambda_output )

  lambda_output = stiff_lambda ( 33.0 )
  print ( "" )
  print ( "  Calling with input and output returns the new default:" )
  print ( "    lambda_output = stiff_lambda ( 33 ) = ", lambda_output )

if ( __name__ == '__main__' ):
  import matplotlib.pyplot as plt
  import numpy as np

  stiff_lambda_test ( )

  LAMBDA = 4.0
  stiff_lambda ( LAMBDA )
  x = np.linspace ( 0.0, 20.0, 101 )
  y = stiff_solution ( x )
  plt.plot ( x, y, 'b-' )
  plt.grid ( True )
  plt.title ( 'y = stiff_solution(x), LAMBDA = ' + str ( LAMBDA ) )
  plt.savefig ( 'stiff_ode.jpg' )
  plt.show ( )
  plt.close ( )