### Test Sundials.jl Installation

Source: https://github.com/sciml/sundials.jl/blob/master/README.md

This code snippet shows how to run the test suite for the Sundials.jl package using the Pkg package manager. This is a useful step to verify that the package has been installed correctly and is functioning as expected.

```julia
using Pkg
Pkg.test("Sundials")
```

--------------------------------

### Install Sundials.jl Package

Source: https://github.com/sciml/sundials.jl/blob/master/README.md

This code snippet demonstrates how to add the Sundials.jl package to a Julia environment using the Pkg package manager. This is the standard method for installing Julia packages.

```julia
using Pkg
Pkg.add("Sundials")
```

--------------------------------

### Using PresetTimeCallback to Modify Solution During Integration

Source: https://context7.com/sciml/sundials.jl/llms.txt

This snippet demonstrates the use of `PresetTimeCallback` to modify the solution at specific times during the integration. The example enforces positivity on the first state variable (`integrator.u[1] = max(integrator.u[1], 0.0)`) at a preset time point.

```julia
using Sundials
using DiffEqCallbacks

# Assume stiff_problem, initial conditions, and tspan are defined
prob = ODEProblem(stiff_problem, [1.0, 0.0], (0.0, 1.0))

# Modify solution during integration
function affect_modify!(integrator)
    integrator.u[1] = max(integrator.u[1], 0.0)  # Enforce positivity
end
cb = PresetTimeCallback([0.5], affect_modify!)
sol = solve(prob, CVODE_BDF(), callback=cb)
```

--------------------------------

### Solve Lorenz Attractor ODE with CVODE_Adams

Source: https://github.com/sciml/sundials.jl/blob/master/README.md

This example demonstrates solving the Lorenz attractor, a classic chaotic system, using the CVODE_Adams method from the Sundials.jl package. It defines the ODE function, initial conditions, and time span, then solves the ODE problem and plots the results.

```julia
using Sundials
function lorenz(du,u,p,t)
 du[1] = 10.0(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob,CVODE_Adams())
using Plots; plot(sol,vars=(1,2,3))
```

--------------------------------

### Using DiscreteCallback to Terminate Integration in Sundials.jl

Source: https://context7.com/sciml/sundials.jl/llms.txt

Shows how to use a `DiscreteCallback` to stop the integration process when a specific condition is met. The example defines a condition based on the state of the solution (`u[1] < 0.1`) and an `affect!` function that calls `terminate!` to end the simulation.

```julia
using Sundials
using DiffEqCallbacks

# Assume stiff_problem, initial conditions, and tspan are defined
prob = ODEProblem(stiff_problem, [1.0, 0.0], (0.0, 1.0))

# Stop when condition met
condition(u, t, integrator) = u[1] < 0.1
affect!(integrator) = terminate!(integrator)
cb = DiscreteCallback(condition, affect!)
sol = solve(prob, CVODE_BDF(), callback=cb)
```

--------------------------------

### DAE Initialization Algorithms (Julia)

Source: https://context7.com/sciml/sundials.jl/llms.txt

Demonstrates using different initialization algorithms in Sundials.jl's IDA solver to compute consistent initial conditions for Differential-Algebraic Equations (DAEs). Covers `BrownFullBasicInit` for index-1 DAEs and `ShampineCollocationInit` which can work without explicit `differential_vars`.

```julia
using Sundials
using DiffEqBase: BrownFullBasicInit, ShampineCollocationInit, CheckInit

# Define the DAE function
function f(res, du, u, p, t)
    res[1] = du[1] - u[2]
    res[2] = du[2] + u[1]
    res[3] = u[1]^2 + u[2]^2 - 1.0  # Algebraic constraint
end

# Initial conditions that do not satisfy the constraint
u₀ = [0.8, 0.5]
du₀ = [0.0, 0.0]
tspan = (0.0, 10.0)
differential_vars = [true, true, false]

# Using BrownFullBasicInit
prob = DAEProblem(f, du₀, u₀, tspan, differential_vars=differential_vars)
sol = solve(prob, IDA(), initializealg=BrownFullBasicInit())

# Using ShampineCollocationInit (does not require differential_vars)
prob = DAEProblem(f, du₀, u₀, tspan)
sol = solve(prob, IDA(), initializealg=ShampineCollocationInit())
```

--------------------------------

### DAE Initialization Strategies in Sundials.jl

Source: https://context7.com/sciml/sundials.jl/llms.txt

Demonstrates different initialization strategies for DAE problems using the IDA solver. Includes checking initial conditions for consistency, default intelligent initialization, and explicitly controlling internal IDA initialization parameters. It also shows how to bypass initialization by providing consistent initial conditions.

```julia
using Sundials

# Assume f, du₀, u₀, tspan, differential_vars are defined

# CheckInit - only validates, does not modify initial conditions
prob = DAEProblem(f, du₀, u₀, tspan, differential_vars=differential_vars)
try
    sol = solve(prob, IDA(), initializealg=CheckInit())
    @error "Should have failed - initial conditions inconsistent"
catch e
    @show e  # InitialConditionError thrown
end

# DefaultInit - intelligently chooses based on problem
prob = DAEProblem(f, du₀, u₀, tspan, differential_vars=differential_vars)
sol = solve(prob, IDA())  # Uses DefaultInit by default

# Control IDA's internal initialization parameters
sol = solve(prob, IDA(
    nonlinear_convergence_coefficient_ic=0.0033,
    max_num_steps_ic=5,
    max_num_jacs_ic=4,
    max_num_iters_ic=10,
    max_num_backs_ic=100,
    use_linesearch_ic=true,
    init_all=false  # false: correct algebraic vars only, true: correct all
), initializealg=BrownFullBasicInit())

# Provide consistent initial conditions to skip initialization
u₀_consistent = [1.0, 0.0, 0.0]
du₀_consistent = [0.0, 1.0, 0.0]
prob = DAEProblem(f, du₀_consistent, u₀_consistent, tspan)
sol = solve(prob, IDA())  # No initialization needed
```

--------------------------------

### IDA Solvers with Automatic Initial Condition Calculation (Julia)

Source: https://context7.com/sciml/sundials.jl/llms.txt

Demonstrates solving DAEs using `idasol` in Sundials.jl. It shows how to automatically calculate consistent initial conditions when provided values do not satisfy the system's constraints, by specifying `diffstates`. It also covers setting solver tolerances and passing user-defined parameters to the residual function.

```julia
using Sundials

# Define the DAE function
function f(t, y, yp, r, userdata)
    k1, k2, k3 = userdata
    r[1] = -k1*y[1] + k2*y[2]*y[3] - yp[1]
    r[2] = k1*y[1] - k3*y[2]^2 - k2*y[2]*y[3] - yp[2]
    r[3] = y[1] + y[2] + y[3] - 1.0
end

# Problem setup
# y0, yp0, t are assumed to be defined elsewhere
differential_vars = [true, true, false]

# Solve with automatic initial condition calculation
yout, ypout = idasol(f, y0, yp0, t, diffstates=differential_vars)

# Solve with consistent initial conditions (no initialization needed)
y0_consistent = [1.0, 0.0, 0.0]
yp0_consistent = [-0.04, 0.04, 0.0]
yout, ypout = idasol(f, y0_consistent, yp0_consistent, t)

# Configure solver tolerances
yout, ypout = idasol(f, y0, yp0, t,
                      reltol=1e-6,
                      abstol=1e-8,
                      diffstates=differential_vars)

# Pass user data to residual function
params = (0.04, 1e4, 3e7)
yout, ypout = idasol(f_with_params, y0, yp0, t, params,
                      diffstates=differential_vars)
```

--------------------------------

### Configuring ARKODE Solver Parameters in Sundials.jl

Source: https://context7.com/sciml/sundials.jl/llms.txt

Illustrates the configuration of the ARKODE solver for ODE problems, focusing on adaptive step size control and predictor methods. Parameters include the order of the RK method, dense output order, predictor method type, convergence coefficients, adaptivity methods, and step size adjustment factors.

```julia
using Sundials

# Assume stiff_problem, initial conditions, and tspan are defined
prob = ODEProblem(stiff_problem, [1.0, 0.0], (0.0, 1.0))

# Configure ARKODE adaptive parameters
sol = solve(prob, ARKODE(
    order=5,                          # RK method order
    dense_order=3,                    # Interpolation polynomial order
    predictor_method=0,               # 0=trivial, 1-4=various predictors
    nonlinear_convergence_coefficient=0.1,  # Convergence threshold
    adaptivity_method=0,              # Adaptivity controller type
    set_optimal_params=false,         # Use default vs optimal stability params
    crdown=0.3,                       # Step decrease factor after failure
    dgmax=0.2,                        # Growth bound on first step
    rdiv=2.3,                         # Divergence threshold
    msbp=20                           # Max steps between predictor evals
))
```

--------------------------------

### Linear Solver Selection for Implicit Methods (Julia)

Source: https://context7.com/sciml/sundials.jl/llms.txt

Illustrates how to select and configure different linear solvers within the CVODE_BDF solver for implicit ODE/DAE integration in Sundials.jl. This includes dense, banded, iterative (GMRES, BCG, PCG, TFQMR), and sparse direct (KLU) solvers, along with options for preconditioners.

```julia
using Sundials
using SparseArrays

# Define ODE problem (example)
prob = ODEProblem((du,u,p,t) -> du .= -u, ones(50), (0.0, 1.0))

# Dense solver
sol = solve(prob, CVODE_BDF(linear_solver=:Dense))

# Banded solver (e.g., tridiagonal Jacobian)
function f_banded(du, u, p, t)
    n = length(u)
    du[1] = -2u[1] + u[2]
    for i in 2:n-1
        du[i] = u[i-1] - 2u[i] + u[i+1]
    end
    du[n] = u[n-1] - 2u[n]
end
prob = ODEProblem(f_banded, ones(100), (0.0, 1.0))
sol = solve(prob, CVODE_BDF(linear_solver=:Band, jac_upper=1, jac_lower=1))

# LapackDense solver (multithreaded BLAS)
prob = ODEProblem((du,u,p,t) -> du .= -u, ones(500), (0.0, 1.0))
sol = solve(prob, CVODE_BDF(linear_solver=:LapackDense))

# GMRES iterative solver
prob = ODEProblem((du,u,p,t) -> du .= -u, ones(10000), (0.0, 1.0))
sol = solve(prob, CVODE_BDF(linear_solver=:GMRES, krylov_dim=30))

# KLU sparse direct solver
function f_sparse(du, u, p, t)
    n = length(u)
    du[1] = -2u[1] + u[2]
    for i in 2:n-1
        du[i] = u[i-1] - 2u[i] + u[i+1]
    end
    du[n] = u[n-1] - 2u[n]
end
jac_prototype = sparse(Tridiagonal(ones(99), -2ones(100), ones(99)))
prob = ODEProblem(ODEFunction(f_sparse, jac_prototype=jac_prototype),
                  ones(100), (0.0, 1.0))
sol = solve(prob, CVODE_BDF(linear_solver=:KLU))

# GMRES with preconditioner
function prec(z, r, p, t, y, fy, gamma, delta, lr)
    z .= r ./ (1.0 .+ gamma)
end
function psetup(p, t, u, du, jok, jcurPtr, gamma)
    jcurPtr[] = true
end
sol = solve(prob, CVODE_BDF(linear_solver=:GMRES,
                            prec=prec,
                            psetup=psetup,
                            prec_side=1))

# Other iterative solvers
sol = solve(prob, CVODE_BDF(linear_solver=:FGMRES, krylov_dim=50))
sol = solve(prob, CVODE_BDF(linear_solver=:BCG))
sol = solve(prob, CVODE_BDF(linear_solver=:PCG))
sol = solve(prob, CVODE_BDF(linear_solver=:TFQMR))
```

--------------------------------

### KINSOL Nonlinear Solver API

Source: https://github.com/sciml/sundials.jl/blob/master/README.md

This documentation describes the high-level API for the KINSOL nonlinear solver within Sundials.jl. It specifies the function signature, arguments, and the expected behavior of the user-provided residual function `f(res,y)`. It also mentions recommended linear solvers and advises on using NLsolve.jl for more robust solutions.

```julia
kinsol(f, y0::Vector{Float64}; 
                userdata::Any = nothing,
                linear_solver=:Dense, jac_upper=0, jac_lower=0)
```

--------------------------------

### Solve Nonlinear Systems with KINSOL in Sundials.jl

Source: https://context7.com/sciml/sundials.jl/llms.txt

Demonstrates solving nonlinear algebraic systems using the KINSOL interface in Sundials.jl. It covers various options including default solvers, line search globalization, iterative solvers like GMRES, banded Jacobians, and tolerance configuration.

```julia
using Sundials

# Simple nonlinear system: solve for roots of system
function f(fy, y)
    fy[1] = y[1]^2 + y[2]^2 - 1.0  # x² + y² = 1 (circle)
    fy[2] = y[2] - y[1]^2           # y = x² (parabola)
end

# Initial guess
y0 = [0.5, 0.5]

# Solve with default dense solver
y_sol = kinsol(f, y0)
@show y_sol
@show f(zeros(2), y_sol)  # Verify solution

# Get full return information (solution and flag)
y_sol, flag = Sundials.___kinsol(f, y0)
@show flag  # Check convergence status

# Use with line search globalization for robustness
y_sol = kinsol(f, y0, strategy=:LineSearch)

# Use iterative GMRES solver for large systems
y_sol = kinsol(f, y0, linear_solver=:GMRES, krylov_dim=50)

# Use banded Jacobian for structured problems
function f_banded(fy, y)
    n = length(y)
    for i in 1:n
        fy[i] = 2.0*y[i]
        if i > 1
            fy[i] -= y[i-1]
        end
        if i < n
            fy[i] -= y[i+1]
        end
    end
end

y0_large = ones(100)
y_sol = kinsol(f_banded, y0_large,
               linear_solver=:Band,
               jac_upper=1,
               jac_lower=1)

# Configure solver tolerances and limits
y_sol = kinsol(f, y0,
               abstol=1e-8,
               maxiters=100,
               maxsetupcalls=10)
```

--------------------------------

### Solve DAE with Sundials.jl

Source: https://context7.com/sciml/sundials.jl/llms.txt

Demonstrates solving a Differential Algebraic Equation (DAE) problem using the DAEProblem and solve functions from Sundials.jl. It shows automatic initialization for inconsistent initial conditions and iterative solvers.

```julia
function f(res, du, u, p, t)
    res[1] = -0.04u[1] + 1e4 * u[2] * u[3] - du[1]
    res[2] = +0.04u[1] - 3e7 * u[2]^2 - 1e4 * u[2] * u[3] - du[2]
    res[3] = u[1] + u[2] + u[3] - 1.0  # Algebraic constraint
end

u₀ = [1.0, 0.0, 0.0]
du₀ = [-0.04, 0.04, 0.0]
tspan = (0.0, 1e5)

# Specify which variables are differential (true) vs algebraic (false)
differential_vars = [true, true, false]
prob = DAEProblem(f, du₀, u₀, tspan, differential_vars=differential_vars)

# Solve with automatic initialization for inconsistent initial conditions
sol = solve(prob, IDA(), initializealg=BrownFullBasicInit())
@show sol.retcode  # Check if solve succeeded

# Provide consistent initial conditions (skip initialization)
u₀_consistent = [1.0, 0.0, 0.0]
du₀_consistent = [-0.04, 0.04, 0.0]
prob_consistent = DAEProblem(f, du₀_consistent, u₀_consistent, tspan,
                             differential_vars=differential_vars)
sol = solve(prob_consistent, IDA())

# Use iterative solver with preconditioner for large DAEs
function prec(z, r, p, t, y, fy, resid, gamma, delta)
    # Simple diagonal preconditioner
    z .= r ./ (1.0 .+ gamma)
end

function psetup(p, t, resid, u, du, gamma)
    # Setup preconditioner (called periodically)
    return nothing
end

sol_iter = solve(prob, IDA(linear_solver=:GMRES, prec=prec, psetup=psetup))

# Configure DAE-specific tolerances and limits
sol = solve(prob, IDA(max_order=3,
                      max_nonlinear_iters=4,
                      nonlinear_convergence_coefficient=0.1),
            reltol=1e-6, abstol=1e-8)
```

--------------------------------

### Enabling Verbose Output for Debugging in Sundials.jl

Source: https://context7.com/sciml/sundials.jl/llms.txt

Shows how to enable verbose output during the solving process for debugging purposes. Setting `verbose=true` in the `solve` function will print step-by-step information, which can be helpful for understanding solver behavior and diagnosing issues.

```julia
using Sundials

# Assume stiff_problem, initial conditions, and tspan are defined
prob = ODEProblem(stiff_problem, [1.0, 0.0], (0.0, 1.0))

# Configure verbosity for debugging
sol = solve(prob, CVODE_BDF(), verbose=true)  # Print step information
```

--------------------------------

### Direct CVODE Interface for ODE Integration in Sundials.jl

Source: https://context7.com/sciml/sundials.jl/llms.txt

Provides a low-level interface to the CVODE solver for ordinary differential equations (ODEs) in Sundials.jl. Allows fine-grained control over solver configuration, including integrator choice, callbacks, user data, and tolerances.

```julia
using Sundials

# Low-level cvode function for time-series integration
function f(t, y, yp)
    yp[1] = -0.5 * y[1]
    yp[2] = -2.0 * y[2]
end

y0 = [1.0, 1.0]
t = 0.0:0.1:5.0

# Solve with default BDF + Dense solver
yout = cvode(f, y0, collect(t))
@show size(yout)  # (51, 2) - rows are time points, columns are variables

# Solve with Adams method for non-stiff problems
yout = cvode(f, y0, collect(t), integrator=:Adams)

# Use callback to monitor/control integration
function callback(mem, t, y)
    if y[1] < 0.1
        println("Variable 1 dropped below 0.1 at t=$t")
        return false  # Stop integration
    end
    return true  # Continue
end

yout = cvode(f, y0, collect(t), callback=callback)

# Pass user data to function
mutable struct MyParams
    decay_rate::Float64
    count::Int
end

function f_with_data(t, y, yp, userdata)
    userdata.count += 1
    yp[1] = -userdata.decay_rate * y[1]
end

params = MyParams(0.5, 0)
yout = cvode(f_with_data, y0, collect(t), params)
@show params.count  # Number of function evaluations

# Configure tolerances
yout = cvode(f, y0, collect(t), reltol=1e-8, abstol=1e-10)
```

--------------------------------

### Configuring Time Stepping Behavior in Sundials.jl Solvers

Source: https://context7.com/sciml/sundials.jl/llms.txt

Demonstrates how to control the time stepping behavior of Sundials.jl solvers (specifically CVODE_BDF shown here) using parameters within the `solve` function. This includes setting the initial step size, minimum and maximum allowed step sizes, and the maximum number of iterations.

```julia
using Sundials

# Assume stiff_problem, initial conditions, and tspan are defined
prob = ODEProblem(stiff_problem, [1.0, 0.0], (0.0, 1.0))

# Configure time stepping behavior through solve parameters
sol = solve(prob, CVODE_BDF(),
    dt=1e-4,                          # Initial step size
    dtmin=1e-10,                      # Minimum allowed step
    dtmax=0.1,                        # Maximum allowed step
    maxiters=1e6,                     # Maximum number of steps
    force_dtmin=true                  # Continue even if dtmin reached
)
```

--------------------------------

### ARKODE: Runge-Kutta Solver for ODEs with Sundials.jl

Source: https://context7.com/sciml/sundials.jl/llms.txt

Utilizes ARKODE for solving ODEs with explicit, implicit, or IMEX (implicit-explicit) Runge-Kutta methods. Supports various orders and specific RK tableaus. For IMEX, requires a SplitODEProblem. Needs Sundials.jl and DifferentialEquations.jl.

```julia
using Sundials
using DifferentialEquations

# Standard ODE with ARKODE
function f(du, u, p, t)
    du[1] = -0.5 * u[1]
    du[2] = -2.0 * u[2]
end

prob = ODEProblem(f, [1.0, 1.0], (0.0, 5.0))

# Use explicit 4th order method
sol = solve(prob, ARKODE(Sundials.Explicit(), order=4))

# Use implicit SDIRK method for stiff problems
sol_implicit = solve(prob, ARKODE(Sundials.Implicit(), order=5))

# IMEX formulation for split ODEs: My' = f_e(t,y) + f_i(t,y)
function f_explicit(du, u, p, t)
    du[1] = u[2]  # Non-stiff advection term
    du[2] = 0.0
end

function f_implicit(du, u, p, t)
    du[1] = 0.0
    du[2] = -100.0 * u[2]  # Stiff diffusion term
end

# Create SplitODEProblem for IMEX
prob_split = SplitODEProblem(f_explicit, f_implicit, [1.0, 1.0], (0.0, 1.0))
sol_imex = solve(prob_split, ARKODE())

# Use specific Runge-Kutta tableau
sol_tableau = solve(prob, ARKODE(Sundials.Explicit(),
                    etable=Sundials.DORMAND_PRINCE_7_4_5))
```

--------------------------------

### Configuring CVODE_BDF Solver Parameters in Sundials.jl

Source: https://context7.com/sciml/sundials.jl/llms.txt

This snippet shows how to fine-tune the CVODE_BDF solver for stiff ODE problems. It demonstrates setting parameters like maximum order, error test failures, nonlinear iterations, convergence failures, and warnings for roundoff errors in time stepping.

```julia
using Sundials

# Assume stiff_problem, initial conditions, and tspan are defined
prob = ODEProblem(stiff_problem, [1.0, 0.0], (0.0, 1.0))

# Configure convergence parameters
sol = solve(prob, CVODE_BDF(
    max_order=3,                      # Limit BDF order (default 5)
    max_error_test_failures=10,       # Max error test fails before reducing step
    max_nonlinear_iters=5,            # Max Newton iterations per step
    max_convergence_failures=15,      # Max convergence fails before stopping
    max_hnil_warns=5                  # Warnings for t+h=t roundoff
))
```

--------------------------------

### Direct IDA Interface for DAE Problems in Sundials.jl

Source: https://context7.com/sciml/sundials.jl/llms.txt

Presents the low-level IDA interface in Sundials.jl for solving Differential Algebraic Equations (DAEs). This interface allows direct access to IDA solver capabilities, including residual formulation and automatic consistent initial condition computation.

```julia
using Sundials

# DAE residual function: F(t, y, yp, r) where F = r
function f(t, y, yp, r)
    r[1] = -0.04y[1] + 1e4*y[2]*y[3] - yp[1]
    r[2] = 0.04y[1] - 3e7*y[2]^2 - 1e4*y[2]*y[3] - yp[2]
    r[3] = y[1] + y[2] + y[3] - 1.0
end

y0 = [1.0, 0.0, 0.0]
yp0 = [-0.04, 0.04, 0.0]
t = [0.0, 1e-4, 1e-3, 1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3, 1e4, 1e5]
```

--------------------------------

### CVODE_BDF: Solve Stiff ODEs with Sundials.jl

Source: https://context7.com/sciml/sundials.jl/llms.txt

Solves stiff ordinary differential equations using the CVODE_BDF solver from Sundials.jl. It supports various linear solvers (Dense, GMRES, Band) and allows custom tolerances. Requires the Sundials.jl and DifferentialEquations.jl packages.

```julia
using Sundials
using DifferentialEquations

# Define the Lorenz system (stiff chaotic ODE)
function lorenz(du, u, p, t)
    du[1] = 10.0(u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8/3) * u[3]
end

u0 = [1.0, 0.0, 0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz, u0, tspan)

# Solve with default Newton + Dense solver
sol = solve(prob, CVODE_BDF())

# Access solution at specific times
@show sol(50.0)  # Interpolated solution at t=50
@show sol.t      # Time points
@show sol.u      # Solution values

# Solve with iterative GMRES solver for large systems
sol_gmres = solve(prob, CVODE_BDF(linear_solver=:GMRES, krylov_dim=30))

# Solve with banded Jacobian (for problems with local coupling)
sol_band = solve(prob, CVODE_BDF(linear_solver=:Band, jac_upper=3, jac_lower=3))

# Solve with custom tolerances
sol_precise = solve(prob, CVODE_BDF(), reltol=1e-8, abstol=1e-10)
```

--------------------------------

### Enabling Stability Limit Detection in CVODE_BDF

Source: https://context7.com/sciml/sundials.jl/llms.txt

This code snippet shows how to enable the stability limit detection feature within the CVODE_BDF solver. This is particularly useful for stiff problems where numerical stability can be an issue, helping the solver to automatically adjust parameters to maintain stability.

```julia
using Sundials

# Assume stiff_problem, initial conditions, and tspan are defined
prob = ODEProblem(stiff_problem, [1.0, 0.0], (0.0, 1.0))

# Enable stability limit detection for BDF
sol = solve(prob, CVODE_BDF(stability_limit_detect=true))
```

--------------------------------

### CVODE_Adams: Solve Non-stiff ODEs with Sundials.jl

Source: https://context7.com/sciml/sundials.jl/llms.txt

Solves non-stiff ordinary differential equations using the CVODE_Adams solver from Sundials.jl. It defaults to functional iteration and supports saving solutions at specified intervals or exact time points. Requires Sundials.jl and DifferentialEquations.jl.

```julia
using Sundials
using DifferentialEquations

# Define a non-stiff problem: simple harmonic oscillator
function harmonic(du, u, p, t)
    ω = 2π
    du[1] = u[2]
    du[2] = -ω^2 * u[1]
end

u0 = [1.0, 0.0]
tspan = (0.0, 10.0)
prob = ODEProblem(harmonic, u0, tspan)

# Solve with Adams method (functional iteration, no linear solver)
sol = solve(prob, CVODE_Adams())

# Save solution at specific intervals
sol = solve(prob, CVODE_Adams(), saveat=0.1)
@show length(sol.t)  # Number of saved points

# Use Newton method instead of functional iteration
sol_newton = solve(prob, CVODE_Adams(method=:Newton, linear_solver=:Dense))

# Specify time stops where solution must be computed exactly
sol = solve(prob, CVODE_Adams(), tstops=[2.5, 5.0, 7.5])
@test all(t in sol.t for t in [2.5, 5.0, 7.5])
```

--------------------------------

### IDA: DAE Solver with Sundials.jl

Source: https://context7.com/sciml/sundials.jl/llms.txt

Solves implicit differential-algebraic equations (DAEs) of the form F(t, y, y') = 0 using the IDA solver from Sundials.jl. It supports index-1 DAEs and automatic or manual specification of variable types. Requires Sundials.jl and DiffEqBase.

```julia
using Sundials
using DiffEqBase: BrownFullBasicInit
```

=== COMPLETE CONTENT === This response contains all available snippets from this library. No additional content exists. Do not make further requests.