USTChEX Documentation - ODE Solvers

USTCHEX_EULER: Explicit (Forward) Euler method for solving first-order ODEs and systems of ODEs

• Description:
  Implements the explicit Euler time-marching method to numerically solve a single first-order ordinary differential equation or a system of coupled first-order ODEs defined directly in Excel cells. The function parses the ODE expressions from worksheet formulas, replaces symbolic variables with their corresponding cell references, and advances the solution from an initial point to a final independent-variable value using a fixed number of uniform steps. The output is a table containing the independent variable in the first column and the computed dependent variable(s) in subsequent columns.

• Syntax:
  =USTCHEX_EULER(Rng, xo, xf, yo, n)

• Rng:
  Range containing the ODE formula(s).
  • For a single ODE, this is one cell containing a formula of the form =f(x,y)
  • For a system of ODEs, this is a vertical range where each cell contains one equation
  Each cell must contain a formula of the form =f(x,y₁,y₂,…). The equal sign is removed internally.

• xo:
  Cell containing the initial value of the independent variable (e.g., time or spatial coordinate). The cell reference or defined name is used symbolically inside the ODE expressions.

• xf:
  Final value of the independent variable at which the integration stops.

• yo:
  Cell or vertical range containing the initial condition(s) of the dependent variable(s).
  • Single cell for a single ODE
  • Vertical range for a system of ODEs
  Cell references or defined names are used symbolically in the ODE expressions.

• n:
  Number of Euler steps to be taken between xo and xf. The step size is computed internally as
  h = (xf − xo) / n

• Returns:
  A 2-D array suitable for spilling into a worksheet.
  • Column 1: independent variable values
  • Columns 2 onward: dependent variable values (one column per ODE)
  The first row corresponds to the initial condition.

• Notes:
  • Uses a fixed step size and the explicit Euler update yₖ₊₁ = yₖ + h·f(xₖ, yₖ)
  • No stability checking or adaptive step sizing is performed
  • Accuracy is first-order in step size
  • Designed for instructional, exploratory, and low-stiffness problems
  • Supports symbolic variable names via defined names or direct cell references

• Example:
  =USTCHEX_EULER(A1:A2, B1, 10, C1:C2, 100)

USTCHEX_BKDEULER: Implicit (Backward) Euler method using Newton–Raphson iteration

• Description:
  Solves a single first-order ordinary differential equation or a system of coupled first-order ODEs using the implicit (Backward Euler) time integration scheme. At each time step, the resulting nonlinear algebraic equation(s) are solved using a Newton–Raphson method with numerically approximated Jacobians. The function supports both scalar and vector systems, optional central or forward finite-difference Jacobian evaluation, adaptive damping for Newton updates, and user-controlled convergence parameters. Output is a tabulated solution history suitable for direct Excel spill.

• Syntax:
  =USTCHEX_BKDEULER(Rng, xo, xf, yo, n, [ctr_chk], [beta_NR], [h_NR], [n_NR], [tol_NR])

• Rng:
  Range containing the ODE formula(s).
  • For a single ODE, this is one cell containing a formula of the form =f(x,y)
  • For a system of ODEs, this is a vertical range where each cell contains one equation
  Each cell must contain a formula of the form =f(x,y₁,y₂,…). The equal sign is removed internally.

• xo:
  Cell containing the initial value of the independent variable (e.g., time). Its cell reference or defined name is symbolically substituted into the ODE expressions.

• xf:
  Final value of the independent variable.

• yo:
  Cell or vertical range containing the initial condition(s) of the dependent variable(s).
  • Single cell for one ODE
  • Vertical range for a system
  Defined names or cell references are both supported.

• n:
  Number of integration steps. The implicit Euler step size is computed as
  h = (xf − xo) / n

• ctr_chk (Optional):
  Logical flag controlling Jacobian evaluation method.
  • TRUE → central finite difference
  • FALSE or omitted → forward finite difference
  Central differencing improves accuracy at the cost of additional evaluations.

• beta_NR (Optional):
  Damping factor applied to Newton–Raphson updates (used for convergence control).
  Default: 1 

• h_NR (Optional):
  Perturbation size used for numerical differentiation of the Jacobian.
  Default: 1×10⁻¹²

• n_NR (Optional):
  Maximum number of Newton–Raphson iterations per time step.
  Default: 512

• tol_NR (Optional):
  Convergence tolerance based on the squared residual norm or update magnitude.
  Default: 1×10⁻¹²

• Returns:
  A 2-D array containing the numerical solution.
  • Column 1: independent variable values
  • Columns 2 onward: dependent variable values
  The first row corresponds to the initial condition.

• Notes:
  • Uses fully implicit Backward Euler discretization
  • Nonlinear systems are solved via Newton–Raphson with numerically constructed Jacobians
  • Adaptive damping improves robustness for stiff or strongly nonlinear systems
  • Suitable for stiff ODEs where explicit Euler is unstable
  • No adaptive time stepping; stability is controlled via implicitness and Newton convergence
  • Designed for research and exploratory numerical analysis in Excel

• Example:
  =USTCHEX_BKDEULER(A1:A3, B1, 20, C1:C3, 200, TRUE)

USTCHEX_RK4: Fixed-step classical fourth-order Runge–Kutta (RK4) solver

• Description:
  Solves a single first-order ordinary differential equation or a system of coupled first-order ODEs using the classical fourth-order explicit Runge–Kutta (RK4) method. The solver advances the solution using a fixed step size computed from the specified interval and number of steps, evaluating four stage derivatives (k₁–k₄) per step. The method supports both scalar equations and vector systems defined directly in Excel cells and returns the complete solution trajectory as a spillable table.

• Syntax:
  =USTCHEX_RK4(Rng, xo, xf, yo, n)

• Rng:
  Range containing the ODE formula(s).
  • For a single ODE, this is one cell containing a formula of the form =f(x,y)
  • For a system of ODEs, this is a vertical range where each cell contains one equation
  Each cell must contain a formula of the form =f(x,y₁,y₂,…). The equal sign is removed internally.

• xo:
  Cell containing the initial value of the independent variable. Its cell reference or defined name is substituted symbolically into the ODE expressions.

• xf:
  Final value of the independent variable.

• yo:
  Cell or vertical range containing the initial condition(s) of the dependent variable(s).

• n:
  Number of fixed integration steps between xo and xf. The step size is computed internally as
  h = (xf − xo) / n

• Returns:
  A 2-D array containing the numerical solution history.
  • Column 1: independent variable values
  • Columns 2 onward: dependent variable values
  The first row corresponds to the initial condition.

• Notes:
  • Explicit fourth-order Runge–Kutta method with fixed step size
  • No adaptive step-size control or error estimation
  • Requires four function evaluations per step
  • Efficient and accurate for smooth, non-stiff problems
  • May become unstable for stiff systems or excessively large step sizes

• Example:
  =USTCHEX_RK4(A1:A3, B1, 10, C1:C3, 200)

USTCHEX_RKF45: Adaptive-step explicit Dormand–Prince version of the Runge–Kutta–Fehlberg (5th-4th order) solver

• Description:
  Solves a single first-order ordinary differential equation or a system of coupled first-order ODEs using the explicit Dormand–Prince embedded Runge–Kutta–Fehlberg 5(4) scheme, the same formulation used by MATLAB’s ode45. At each integration step, both fifth-order and embedded fourth-order solutions are computed from a shared set of stage evaluations, and their difference is used to estimate the local truncation error. The step size is adjusted dynamically to satisfy a user-defined tolerance. The solver supports scalar and vector systems, enforces minimum and maximum step-size limits, and applies local polynomial interpolation to exactly match the final independent-variable value when overshoot occurs.

• Syntax:
  =USTCHEX_RKF45(Rng, xo, xf, yo, n, [tol], [h], [hmin], [hmax], [lim])

• Rng:
  Range containing the ODE formula(s).
  • For a single ODE, this is one cell containing a formula of the form =f(x,y)
  • For a system of ODEs, this is a vertical range where each cell contains one equation
  Each cell must contain a formula of the form =f(x,y₁,y₂,…). The equal sign is removed internally.

• xo:
  Cell containing the initial value of the independent variable. Its cell reference or defined name is substituted symbolically into the ODE expressions.

• xf:
  Final value of the independent variable.

• yo:
  Cell or vertical range containing the initial condition(s) of the dependent variable(s).

• n:
  Maximum number of solution points to be generated. The solver may terminate earlier if the final independent-variable value is reached.

• tol (Optional):
  Local error tolerance used for adaptive step-size control based on the difference between the fourth- and fifth-order solutions.
  Default: 1×10⁻⁸

• h (Optional):
  Initial step size for the integration procedure.
  Default: 0.1

• hmin (Optional):
  Minimum allowable step size. Integration halts if further reduction is required beyond this limit.
  Default: 1×10⁻¹⁵

• hmax (Optional):
  Maximum allowable step size.
  Default: 1

• lim (Optional):
  Maximum number of consecutive step-size reductions permitted before termination.
  Default: 1000

• Returns:
  A 2-D array containing the numerical solution history.
  • Column 1: independent variable values
  • Columns 2 onward: dependent variable values
  The solution is interpolated to end exactly at xf if the adaptive step overshoots the final value.

• Notes:
  • Explicit Runge–Kutta–Fehlberg 4(5) method with embedded error estimation
  • Automatic step-size adaptation improves efficiency for smooth, non-stiff problems
  • Requires six to seven function evaluations per accepted step
  • Not suitable for stiff systems due to explicit formulation
  • Enforces user-defined minimum and maximum step sizes
  • Uses local polynomial interpolation to exactly match the final independent-variable value

• Example:
  =USTCHEX_RKF45(A1:A3, B1, 20, C1:C3, 500)

USTCHEX_BDF2: Second-order Backward Differentiation Formula (BDF2) with implicit Newton–Raphson solution\

• Description:
  Solves a single first-order ordinary differential equation or a system of coupled first-order ODEs using the second-order Backward Differentiation Formula (BDF2). The method is fully implicit and requires the solution of nonlinear algebraic equations at each time step, which is accomplished via Newton–Raphson iteration with numerically approximated Jacobians and adaptive damping. A Backward Euler startup phase is used internally to generate the initial history required by the multistep scheme. The function supports scalar and vector systems and returns the full solution trajectory as a spillable table.

• Syntax:
  =USTCHEX_BDF2(Rng, xo, xf, yo, n_BDF2, [n_EULER], [ctr_chk], [beta_NR], [h_NR], [n_NR], [tol_NR])

• Rng:
  Range containing the ODE formula(s).
  • For a single ODE, this is one cell containing a formula of the form =f(x,y)
  • For a system of ODEs, this is a vertical range where each cell contains one equation
  Each cell must contain a formula of the form =f(x,y₁,y₂,…). The equal sign is removed internally.

• xo:
  Cell containing the initial value of the independent variable. Its cell reference or defined name is substituted symbolically into the ODE expressions.

• xf:
  Final value of the independent variable.

• yo:
  Cell or vertical range containing the initial condition(s) of the dependent variable(s).

• n_BDF2:
  Number of BDF2 time steps between xo and xf. The BDF2 step size is computed internally as
  h_BDF2 = (xf − xo) / n_BDF2

• n_EULER (Optional):
  Number of Backward Euler substeps used to generate the initial solution history for BDF2.
  Default: 10

• ctr_chk (Optional):
  Logical flag controlling Jacobian evaluation.
  • TRUE → central finite differences
  • FALSE or omitted → forward finite differences

• beta_NR (Optional):
  Damping factor applied to Newton–Raphson updates (used for convergence control).
  Default: 1 

• h_NR (Optional):
  Perturbation size used for numerical Jacobian evaluation.
  Default: 0.01

• n_NR (Optional):
  Maximum number of Newton–Raphson iterations per time step.
  Default: 512

• tol_NR (Optional):
  Convergence tolerance based on the squared residual norm.
  Default: 1×10⁻¹²

• Returns:
  A 2-D array containing the numerical solution history.
  • Column 1: independent variable values
  • Columns 2 onward: dependent variable values
  The first row corresponds to the initial condition.

• Notes:
  • Uses second-order implicit BDF time discretization
  • Backward Euler startup is automatically performed to initialize the multistep scheme
  • Newton–Raphson iteration employs numerically constructed Jacobians and adaptive damping
  • Well-suited for stiff and moderately stiff systems
  • Fixed time step; stability is achieved through implicit formulation rather than adaptivity

• Example:
  =USTCHEX_BDF2(A1:A3, B1, 50, C1:C3, 200)

USTCHEX_BDF4: Fourth-order Backward Differentiation Formula (BDF4) with multistage implicit startup and Newton–Raphson solution

Description:
Solves a single first-order ordinary differential equation or a system of coupled first-order ODEs using the fourth-order Backward Differentiation Formula (BDF4). Because BDF4 is a multistep implicit method, the solver automatically constructs the required solution history through a hierarchical startup sequence: Backward Euler (BDF1), BDF2, and BDF3, each solved implicitly via Newton–Raphson iteration. Numerical Jacobians are computed using forward or central finite differences, and damped Newton updates are applied to improve nonlinear convergence. The function returns the complete solution trajectory as a spillable table.

• Syntax:
  =USTCHEX_BDF4(Rng, xo, xf, yo, n_BDF4, [n_EULER], [n_BDF2], [n_BDF3], [ctr_chk], [beta_NR], [h_NR], [n_NR], [tol_NR])

• Rng:
  Range containing the ODE formula(s).
  • Single cell for a scalar ODE
  • Vertical range for systems of ODEs
  Each cell must contain a formula of the form =f(x,y₁,y₂,…). The leading equal sign is removed internally.

• xo:
  Cell containing the initial value of the independent variable. The cell reference or defined name is substituted symbolically into the ODE expressions.

• xf:
  Final value of the independent variable.

• yo:
  Cell or vertical range containing the initial condition(s) of the dependent variable(s).

• n_BDF4:
  Number of BDF4 time steps between xo and xf. The BDF4 step size is computed as
  h₄ = (xf − xo) / n_BDF4

• n_EULER (Optional):
  Number of Backward Euler substeps used to generate the initial solution history.
  Default: 10

• n_BDF2 (Optional):
  Number of BDF2 substeps used after the Euler startup.
  Default: 10

• n_BDF3 (Optional):
  Number of BDF3 substeps used after the BDF2 startup.
  Default: 10

• ctr_chk (Optional):
  Logical flag controlling Jacobian evaluation method.
  • TRUE → central finite differences
  • FALSE or omitted → forward finite differences

• beta_NR (Optional):
  Damping factor applied to Newton–Raphson updates (used for convergence control).
  Default: 1

• h_NR (Optional):
  Perturbation size used for numerical Jacobian evaluation.
  Default: 1×10⁻⁴

• n_NR (Optional):
  Maximum number of Newton–Raphson iterations per time step.
  Default: 512

• tol_NR (Optional):
  Convergence tolerance based on the squared residual norm.
  Default: 1×10⁻¹²

• Returns:
  A 2-D array containing the numerical solution history.
  • Column 1: independent variable values
  • Columns 2 onward: dependent variable values
  The first row corresponds to the initial condition.

• Notes:
  • Uses fourth-order implicit BDF time discretization
  • Automatically performs Euler → BDF2 → BDF3 → BDF4 startup sequencing
  • Nonlinear systems are solved using Newton–Raphson with numerically constructed Jacobians
  • Fixed step size throughout; stability derives from implicit formulation and multistep history
  • Suitable for stiff and strongly coupled systems when higher-order accuracy is required

• Example:
  =USTCHEX_BDF4(A1:A3, B1, 100, C1:C3, 400)

USTCHEX_BDF12: Adaptive-step implicit BDF1–BDF2 solver with adaptive step-size control and Newton–Raphson iteration

• Description:
  Solves a single first-order ordinary differential equation or a system of coupled first-order ODEs using an implicit Backward Differentiation framework that combines BDF1 (Backward Euler) and BDF2 with automatic step-size control. At each macro step, a predictor–corrector strategy compares BDF2 and Euler-consistent solutions to estimate the local truncation error and adjust the step size accordingly. The resulting nonlinear systems are solved using Newton–Raphson iteration with numerically approximated Jacobians and damped updates to improve convergence robustness. The solver supports scalar and vector systems, enforces minimum and maximum step-size limits, and performs a hard snap to the final independent-variable value using local polynomial interpolation if overshoot occurs.

• Syntax:
  =USTCHEX_BDF12(Rng, xo, xf, yo, n, [h_ODE], [hmin_ODE], [hmax_ODE], [tol_ODE], [lim_ODE], [ctr_chk], [beta_NR], [h_NR], [n_NR], [tol_NR])

• Rng:
  Range containing the ODE formula(s).
  • For a single ODE, this is one cell containing a formula of the form =f(x,y)
  • For a system of ODEs, this is a vertical range where each cell contains one equation
  Each cell must contain a formula of the form =f(x,y₁,y₂,…). The equal sign is removed internally.

• xo:
  Cell containing the initial value of the independent variable. Its cell reference or defined name is symbolically substituted into the ODE expressions.

• xf:
  Final value of the independent variable.

• yo:
  Cell or vertical range containing the initial condition(s) of the dependent variable(s).

• n:
  Maximum number of accepted solution points to be generated. The solver may terminate earlier if the final independent-variable value is reached.

• h_ODE (Optional):
  Initial step size for the adaptive integration procedure.
  Default: 0.1

• hmin_ODE (Optional):
  Minimum allowable step size. Integration halts if further reduction is required beyond this limit.
  Default: 1×10⁻¹⁵

• hmax_ODE (Optional):
  Maximum allowable step size.
  Default: 100

• tol_ODE (Optional):
  Local error tolerance used for step-size adaptation based on the difference between Euler and BDF2 solutions.
  Default: 1×10⁻⁵

• lim_ODE (Optional):
  Maximum number of consecutive step-size reductions permitted before termination.
  Default: 1000

• ctr_chk (Optional):
  Logical flag controlling Jacobian evaluation method in Newton–Raphson iterations.
  • TRUE → central finite differences
  • FALSE or omitted → forward finite differences

• beta_NR (Optional):
  Damping factor applied to Newton–Raphson updates for convergence control.
  Default: 1

• h_NR (Optional):
  Perturbation size used for numerical Jacobian evaluation.
  Default: 0.01

• n_NR (Optional):
  Maximum number of Newton–Raphson iterations per implicit solve.
  Default: 512

• tol_NR (Optional):
  Convergence tolerance for Newton–Raphson iteration based on the residual norm.
  Default: 1×10⁻¹²

• Returns:
  A 2-D array containing the numerical solution history.
  • Column 1: independent variable values
  • Columns 2 onward: dependent variable values
  The solution is trimmed or interpolated to end exactly at xf unless the iteration count limit n is reached first.

• Notes:
  • Combines BDF1 and BDF2 to enable local error estimation and adaptive step-size control
  • Fully implicit formulation with Newton–Raphson solution at each accepted step
  • Supports stiff and strongly nonlinear systems
  • Enforces user-defined minimum and maximum step sizes
  • Performs controlled interpolation to exactly match the final independent-variable value when overshoot occurs
  • Fixed-order method with variable step size

• Example:
  =USTCHEX_BDF12(A1:A3, B1, 25, C1:C3, 500)

USTCHEX_BDF23: Adaptive-step implicit BDF2–BDF3 solver with local error control and Newton–Raphson iteration

• Description:
  Solves a single first-order ordinary differential equation or a system of coupled first-order ODEs using an implicit variable-step Backward Differentiation framework based on second- and third-order formulas (BDF2–BDF3). At each adaptive macro-step, the solver constructs BDF2 and BDF3 solutions over a short internal history and estimates the local truncation error from their difference. The step size is then adjusted automatically to satisfy a user-specified tolerance. All implicit stages are solved using Newton–Raphson iteration with numerically approximated Jacobians and damped updates to improve nonlinear convergence. The method advances the solution in blocks of three points and includes a hard snap to the final independent-variable value using local polynomial interpolation when overshoot occurs.

• Syntax:
  =USTCHEX_BDF23(Rng, xo, xf, yo, n, [h_ODE], [hmin_ODE], [hmax_ODE], [tol_ODE], [lim_ODE], [ctr_chk], [beta_NR], [h_NR], [n_NR], [tol_NR])

• Rng:
  Range containing the ODE formula(s).
  • For a single ODE, this is one cell containing a formula of the form =f(x,y)
  • For a system of ODEs, this is a vertical range where each cell contains one equation
  Each cell must contain a formula of the form =f(x,y₁,y₂,…). The equal sign is removed internally.

• xo:
  Cell containing the initial value of the independent variable. Its cell reference or defined name is substituted symbolically into the ODE expressions.

• xf:
  Final value of the independent variable.

• yo:
  Cell or vertical range containing the initial condition(s) of the dependent variable(s).

• n:
  Maximum number of solution points to be generated. The solver may terminate earlier if the final independent-variable value is reached.

• h_ODE (Optional):
  Initial step size for the adaptive integration procedure.
  Default: 0.1

• hmin_ODE (Optional):
  Minimum allowable step size. Integration halts if further reduction is required beyond this limit.
  Default: 1×10⁻¹⁵

• hmax_ODE (Optional):
  Maximum allowable step size.
  Default: 100

• tol_ODE (Optional):
  Local error tolerance used for adaptive step-size control based on the difference between BDF3 and BDF2 solutions.
  Default: 1×10⁻⁴

• lim_ODE (Optional):
  Maximum number of consecutive step-size reductions permitted before termination.
  Default: 1000

• ctr_chk (Optional):
  Logical flag controlling Jacobian evaluation method in Newton–Raphson iterations.
  • TRUE → central finite differences
  • FALSE or omitted → forward finite differences

• beta_NR (Optional):
  Damping factor applied to Newton–Raphson updates for convergence control.
  Default: 1

• h_NR (Optional):
  Perturbation size used for numerical Jacobian evaluation.
  Default: 0.01

• n_NR (Optional):
  Maximum number of Newton–Raphson iterations per implicit solve.
  Default: 512

• tol_NR (Optional):
  Convergence tolerance for Newton–Raphson iteration based on the residual norm.
  Default: 1×10⁻¹²

• Returns:
  A 2-D array containing the numerical solution history.
  • Column 1: independent variable values
  • Columns 2 onward: dependent variable values
  The solution is trimmed or interpolated to end exactly at xf unless the iteration count limit n is reached first.

• Notes:
  • Combines BDF2 and BDF3 to enable reliable local error estimation and adaptive step-size control
  • Fully implicit formulation with Newton–Raphson solution at each accepted step
  • Advances the solution in three-point blocks while retaining second-order continuity
  • Enforces user-defined minimum and maximum step sizes
  • Performs controlled interpolation to exactly match the final independent-variable value when overshoot occurs
  • Recommended general-purpose solver for stiff and moderately stiff systems when robustness and efficiency are both required

• Example:
  =USTCHEX_BDF23(A1:A3, B1, 30, C1:C3, 500)

USTCHEX_BDF34: Adaptive-step implicit BDF3–BDF4 solver with local error control and Newton–Raphson iteration

• Description:
  Solves a single first-order ordinary differential equation or a system of coupled first-order ODEs using an implicit variable-step Backward Differentiation framework based on third- and fourth-order formulas (BDF3–BDF4). At each adaptive macro-step, the solver constructs BDF3 and BDF4 solutions over an internally generated multistep history and estimates the local truncation error from their difference. The step size is adjusted automatically to satisfy a user-defined tolerance. All implicit stages are solved using Newton–Raphson iteration with numerically approximated Jacobians and damped updates to improve nonlinear convergence. The method advances the solution in four-point blocks and includes a hard snap to the final independent-variable value using local polynomial interpolation when overshoot occurs.

• Syntax:
  =USTCHEX_BDF34(Rng, xo, xf, yo, n, [h_ODE], [hmin_ODE], [hmax_ODE], [tol_ODE], [lim_ODE], [ctr_chk], [beta_NR], [h_NR], [n_NR], [tol_NR])

• Rng:
  Range containing the ODE formula(s).
  • For a single ODE, this is one cell containing a formula of the form =f(x,y)
  • For a system of ODEs, this is a vertical range where each cell contains one equation
  Each cell must contain a formula of the form =f(x,y₁,y₂,…). The equal sign is removed internally.

• xo:
  Cell containing the initial value of the independent variable. Its cell reference or defined name is substituted symbolically into the ODE expressions.

• xf:
  Final value of the independent variable.

• yo:
  Cell or vertical range containing the initial condition(s) of the dependent variable(s).

• n:
  Maximum number of solution points to be generated. The solver may terminate earlier if the final independent-variable value is reached.

• h_ODE (Optional):
  Initial step size for the adaptive integration procedure.
  Default: 0.1

• hmin_ODE (Optional):
  Minimum allowable step size. Integration halts if further reduction is required beyond this limit.
  Default: 1×10⁻¹⁵

• hmax_ODE (Optional):
  Maximum allowable step size.
  Default: 100

• tol_ODE (Optional):
  Local error tolerance used for adaptive step-size control based on the difference between BDF4 and BDF3 solutions.
  Default: 1×10⁻⁴

• lim_ODE (Optional):
  Maximum number of consecutive step-size reductions permitted before termination.
  Default: 1000

• ctr_chk (Optional):
  Logical flag controlling Jacobian evaluation method in Newton–Raphson iterations.
  • TRUE → central finite differences
  • FALSE or omitted → forward finite differences

• beta_NR (Optional):
  Damping factor applied to Newton–Raphson updates for convergence control.
  Default: 1

• h_NR (Optional):
  Perturbation size used for numerical Jacobian evaluation.
  Default: 0.01

• n_NR (Optional):
  Maximum number of Newton–Raphson iterations per implicit solve.
  Default: 512

• tol_NR (Optional):
  Convergence tolerance for Newton–Raphson iteration based on the residual norm.
  Default: 1×10⁻¹²

• Returns:
  A 2-D array containing the numerical solution history.
  • Column 1: independent variable values
  • Columns 2 onward: dependent variable values
  The solution is trimmed or interpolated to end exactly at xf unless the iteration count limit n is reached first.

• Notes:
  • Combines BDF3 and BDF4 to enable higher-order local error estimation
  • Fully implicit formulation with Newton–Raphson solution at each accepted step
  • Advances the solution in four-point blocks, increasing startup and computational overhead
  • Enforces user-defined minimum and maximum step sizes
  • Performs controlled interpolation to exactly match the final independent-variable value when overshoot occurs
  • Generally less efficient and less robust than BDF23 for most stiff and moderately stiff systems

• Example:
  =USTCHEX_BDF34(A1:A3, B1, 30, C1:C3, 500)