RustedSciThe 0.4.15

//!
//! # BVP_sci_faer - High-Performance Boundary Value Problem Solver
//!
//! This module implements a state-of-the-art 4th order collocation algorithm with residual control,
//! translated and optimized from SciPy's _bvp.py. It provides robust and efficient solutions for
//! boundary value problems in ordinary differential equations using the faer linear algebra library.
//!
//! ## Algorithm Features
//! - **4th Order Collocation**: Uses cubic spline interpolation with C1 continuity
//! - **Adaptive Mesh Refinement**: Automatic grid adaptation based on residual estimates
//! - **Residual Control**: 5-point Lobatto quadrature for accurate error estimation
//! - **Sparse Matrix Operations**: Efficient handling of large, sparse Jacobian matrices
//! - **Newton's Method**: Robust nonlinear solver with backtracking line search
//! - **Parameter Estimation**: Support for unknown parameters in the ODE system
//!
//! ## Performance Optimizations
//! - **faer Linear Algebra**: High-performance sparse matrix operations
//! - **Parallel Jacobian Evaluation**: Concurrent computation of derivatives
//! - **Memory Efficient**: Sparse storage for banded and structured matrices
//! - **Adaptive Strategies**: Smart Jacobian recomputation and step size control
//!
//! ## Supported Problem Types
//! - Linear and nonlinear boundary value problems
//! - Systems of coupled ODEs with up to thousands of equations
//! - Multi-point boundary conditions
//! - Parameter estimation problems
//! - Stiff and non-stiff differential equations
//!
//! ## Function Interface Design
//! Unlike traditional ODE solvers that work with single points, this solver uses vectorized
//! evaluation across the entire mesh for optimal performance:
//!
//! ```rust,ignore
//! // ODE Function: f(x, y, p) -> f_values
//! // x: mesh points (m,)  
//! // y: solution at all mesh points (n×m)
//! // p: parameters (k,)
//! // Returns: RHS values at ALL mesh points (n×m)
//! pub type ODEFunction = dyn Fn(&faer_col, &faer_dense_mat, &faer_col) -> faer_dense_mat;
//!
//! // Jacobian: Returns sparse matrices for each mesh point
//! pub type ODEJacobian = dyn Fn(&faer_col, &faer_dense_mat, &faer_col)
//!                        -> (Vec<faer_mat>, Option<Vec<faer_mat>>);
//! ```
//!
//! ## Main Solver Functions
//! - `solve_bvp()`: Main solver with automatic mesh refinement
//! - `solve_bvp_sparse()`: Optimized version for large sparse systems
//! - `solve_newton()`: Core Newton iteration with line search
//! - `collocation_fun()`: Collocation residual computation
//! - `estimate_rms_residuals()`: Error estimation using Lobatto quadrature
//!
//! ## Usage Example
//! ```rust,ignore
//! // Define ODE: y'' + y = 0
//! let fun = |x: &faer_col, y: &faer_dense_mat, p: &faer_col| {
//!     let mut f = faer_dense_mat::zeros(2, y.ncols());
//!     for j in 0..y.ncols() {
//!         *f.get_mut(0, j) = *y.get(1, j);        // y' = z
//!         *f.get_mut(1, j) = -*y.get(0, j);       // z' = -y  
//!     }
//!     f
//! };
//!
//! // Boundary conditions: y(0) = 0, y(π) = 0
//! let bc = |ya: &faer_col, yb: &faer_col, p: &faer_col| {
//!     faer_col::from_fn(2, |i| match i {
//!         0 => ya[0],     // y(0) = 0
//!         1 => yb[0],     // y(π) = 0
//!         _ => 0.0,
//!     })
//! };
//!
//! let result = solve_bvp(&fun, &bc, x, y_guess, None, None,
//!                        None, None, 1e-6, 1000, 1, None, None)?;
//! ```
//!
//! Boundary value problem solver - Rust translation from SciPy's _bvp.py
//!
//! This module implements a 4th order collocation algorithm with residual control
//! similar to the SciPy BVP solver, translated from Python to Rust.
//! Important notes about construction of functions in this numerical method
//! In most of the others methods Residuals function will be as follows
//! : dyn Fn(f64, &faer_col) -> faer_col - takes argument value and unknowns
//! vector, returns residuals, and Jacobian will be as follows
//! Jacobian: dyn Fn(f64, &faer_col) -> faer_mat - takes argument value and
//!  unknowns vector, returns Jacobian matrix.
//! but here residual fun is arranged as follows:
//! pub type ODEFunction =
//! dyn Fn(&faer_col, &faer_dense_mat, &faer_col) -> faer_dense_mat;
//! Takes (x, y, p) where:
//! x: mesh points (vector)
//! y: solution values at all mesh points (n×m matrix)
//! p: parameters (vector)
//! Returns: RHS values at ALL mesh points (n×m matrix)
//! Jacobian function is arranged as follows
//! pub type ODEJacobian = dyn Fn(&faer_col, &faer_dense_mat, &faer_col) -> (Vec<faer_mat>, Option<Vec<faer_mat>>);
//! Takes same inputs as ODE function
//! Returns: (df_dy, df_dp) where:
//! df_dy: Vector of sparse matrices, one (n×n) matrix per mesh point
//! df_dp: Optional vector of sparse matrices for parameter derivatives
//!

use crate::numerical::BVP_Damp::BVP_utils::CustomTimer;
//use crate::numerical::BVP_sci::BVP_sci_utils::{final_jacobian_diagnostics, size_of_matrix};
use crate::numerical::BVP_sci::BVP_sci_banded::sparse_global_jac_to_banded;
use crate::numerical::BVP_sci::BVP_sci_bordered_banded::{
    BvpSciBandedRoute, BvpSciBandedRoutePolicy, profile_bordered_banded_global_jacobian,
};
use crate::numerical::BVP_sci::BVP_sci_bordered_solver::{
    BvpSciBorderedSolveWorkspace, BvpSciBorderedStructuredFactorization,
    assemble_bordered_banded_blocks_from_pointwise, extract_bordered_banded_blocks,
    factor_bordered_banded_structured,
};
use crate::numerical::optimization::PPoly::{Extrapolate, PPoly};
use crate::somelinalg::banded::LapackStyleBandedLuFaithful;
use crate::somelinalg::banded::banded_assembly::BandedAssembly;
use faer::col::{Col, ColRef};
use faer::linalg::solvers::Solve;
use faer::mat::{Mat, MatRef};
use faer::sparse::{SparseColMat, SymbolicSparseColMat, Triplet};
use log::info;
use nalgebra::{DMatrix, linalg::SVD};

use std::{collections::HashMap, mem::size_of, time::Instant};
pub type faer_mat = SparseColMat<usize, f64>;
pub type faer_col = Col<f64>;
pub type faer_dense_mat = Mat<f64>;

/// Native pointwise banded Jacobian provider.
///
/// The first tuple item contains one state Jacobian per mesh point. Parameter
/// derivatives remain sparse because the current direct bordered route is
/// intentionally limited to parameter-free endpoint problems.
pub type ODEBandedJacobian =
    dyn Fn(&faer_col, &faer_dense_mat, &faer_col) -> (Vec<BandedAssembly>, Option<Vec<faer_mat>>);

#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
pub enum BvpSciLinearSolvePolicy {
    /// Production default: preserve the existing faer sparse LU route.
    #[default]
    Sparse,
    /// Select the native linear route from the global Jacobian structure.
    ///
    /// Compact whole matrices use scalar banded LU. Parameter-free endpoint
    /// systems with a compact collocation body use the bordered-banded solver.
    /// Unsupported layouts and failed native factorizations fall back to Sparse.
    AutoBanded,
    /// Require the whole global Jacobian to be compact enough for scalar banded
    /// LU.  If the route planner rejects it, Newton reports a singular linear
    /// step instead of silently falling back.
    RequireFullBanded,
    /// Require the endpoint-BC bordered-banded route without Sparse fallback.
    ///
    /// This remains useful as a strict diagnostics and acceptance policy even
    /// though `AutoBanded` can now select the same route safely.
    ExperimentalBorderedBanded,
}

#[derive(Clone, Debug, Default, PartialEq, Eq)]
pub struct GlobalJacobianDiagnostics {
    pub nrows: usize,
    pub ncols: usize,
    pub nnz: usize,
    pub dense_equivalent_bytes: usize,
    pub sparse_values_bytes: usize,
    pub sparse_row_index_bytes: usize,
    pub sparse_col_ptr_bytes: usize,
    pub sparse_storage_bytes: usize,
    pub dense_equivalent_kib: usize,
    pub sparse_storage_kib: usize,
    pub dense_equivalent_mib: usize,
    pub sparse_storage_mib: usize,
    pub dense_to_sparse_permille: usize,
    pub nonfinite_values: usize,
    pub structurally_empty_rows: usize,
    pub structurally_empty_cols: usize,
    pub numerically_zero_rows: usize,
    pub numerically_zero_cols: usize,
}

impl GlobalJacobianDiagnostics {
    pub fn is_definitely_singular(&self) -> bool {
        self.nonfinite_values > 0
            || self.structurally_empty_rows > 0
            || self.structurally_empty_cols > 0
            || self.numerically_zero_rows > 0
            || self.numerically_zero_cols > 0
    }
}

pub fn inspect_global_jacobian(mat: &faer_mat) -> GlobalJacobianDiagnostics {
    let (nrows, ncols) = mat.shape();
    let mut row_has_structure = vec![false; nrows];
    let mut col_has_structure = vec![false; ncols];
    let mut row_has_nonzero = vec![false; nrows];
    let mut col_has_nonzero = vec![false; ncols];
    let mut nnz = 0usize;
    let mut nonfinite_values = 0usize;

    for col in 0..ncols {
        let rows = mat.as_dyn().row_idx_of_col(col).collect::<Vec<_>>();
        if !rows.is_empty() {
            col_has_structure[col] = true;
        }
        for row in rows {
            nnz += 1;
            row_has_structure[row] = true;
            let value = mat[(row, col)];
            if !value.is_finite() {
                nonfinite_values += 1;
            }
            if value != 0.0 {
                row_has_nonzero[row] = true;
                col_has_nonzero[col] = true;
            }
        }
    }

    let dense_equivalent_bytes = nrows.saturating_mul(ncols).saturating_mul(size_of::<f64>());
    let sparse_values_bytes = nnz.saturating_mul(size_of::<f64>());
    let sparse_row_index_bytes = nnz.saturating_mul(size_of::<usize>());
    // faer stores sparse matrices in compressed-column form.  This is a stable
    // solver-level estimate: values + row indices + one column pointer per col.
    let sparse_col_ptr_bytes = ncols.saturating_add(1).saturating_mul(size_of::<usize>());
    let sparse_storage_bytes = sparse_values_bytes
        .saturating_add(sparse_row_index_bytes)
        .saturating_add(sparse_col_ptr_bytes);
    let dense_equivalent_kib = bytes_to_kib_ceil(dense_equivalent_bytes);
    let sparse_storage_kib = bytes_to_kib_ceil(sparse_storage_bytes);
    let dense_equivalent_mib = bytes_to_mib_ceil(dense_equivalent_bytes);
    let sparse_storage_mib = bytes_to_mib_ceil(sparse_storage_bytes);
    let dense_to_sparse_permille = if sparse_storage_bytes > 0 {
        dense_equivalent_bytes
            .saturating_mul(1000)
            .saturating_div(sparse_storage_bytes)
    } else {
        0
    };

    GlobalJacobianDiagnostics {
        nrows,
        ncols,
        nnz,
        dense_equivalent_bytes,
        sparse_values_bytes,
        sparse_row_index_bytes,
        sparse_col_ptr_bytes,
        sparse_storage_bytes,
        dense_equivalent_kib,
        sparse_storage_kib,
        dense_equivalent_mib,
        sparse_storage_mib,
        dense_to_sparse_permille,
        nonfinite_values,
        structurally_empty_rows: row_has_structure
            .iter()
            .filter(|has_structure| !**has_structure)
            .count(),
        structurally_empty_cols: col_has_structure
            .iter()
            .filter(|has_structure| !**has_structure)
            .count(),
        numerically_zero_rows: row_has_nonzero
            .iter()
            .filter(|has_nonzero| !**has_nonzero)
            .count(),
        numerically_zero_cols: col_has_nonzero
            .iter()
            .filter(|has_nonzero| !**has_nonzero)
            .count(),
    }
}

fn bytes_to_kib_ceil(bytes: usize) -> usize {
    bytes.saturating_add(1023) / 1024
}

fn bytes_to_mib_ceil(bytes: usize) -> usize {
    bytes.saturating_add(1024 * 1024 - 1) / (1024 * 1024)
}

fn faer_dense_to_dmatrix(mat: &faer_dense_mat) -> DMatrix<f64> {
    DMatrix::from_fn(mat.nrows(), mat.ncols(), |row, col| *mat.get(row, col))
}

fn dmatrix_to_faer_dense(mat: &DMatrix<f64>) -> faer_dense_mat {
    faer_dense_mat::from_fn(mat.nrows(), mat.ncols(), |row, col| mat[(row, col)])
}
#[allow(dead_code)]
fn sparse_to_dmatrix(mat: &faer_mat) -> DMatrix<f64> {
    let (nrows, ncols) = mat.shape();
    let mut dense = DMatrix::zeros(nrows, ncols);
    for col in 0..ncols {
        for row in mat.as_dyn().row_idx_of_col(col) {
            dense[(row, col)] = mat[(row, col)];
        }
    }
    dense
}
#[allow(dead_code)]
fn dmatrix_to_sparse(mat: &DMatrix<f64>) -> faer_mat {
    let mut triplets = Vec::new();
    for col in 0..mat.ncols() {
        for row in 0..mat.nrows() {
            let value = mat[(row, col)];
            if value != 0.0 {
                triplets.push(Triplet::new(row, col, value));
            }
        }
    }
    faer_mat::try_new_from_triplets(mat.nrows(), mat.ncols(), &triplets).unwrap()
}

fn pseudoinverse(mat: &DMatrix<f64>) -> Result<DMatrix<f64>, String> {
    let svd = SVD::new(mat.clone(), true, true);
    let u = svd
        .u
        .ok_or_else(|| "SVD failed to compute left singular vectors".to_string())?;
    let vt = svd
        .v_t
        .ok_or_else(|| "SVD failed to compute right singular vectors".to_string())?;
    let singular_values = svd.singular_values;
    let sigma_max = singular_values
        .iter()
        .fold(0.0_f64, |acc, value| acc.max(*value));
    let tol = (mat.nrows().max(mat.ncols()) as f64) * f64::EPSILON * sigma_max.max(1.0);
    let mut sigma_inv = DMatrix::zeros(vt.nrows(), u.ncols());
    for (idx, sigma) in singular_values.iter().enumerate() {
        if *sigma > tol {
            sigma_inv[(idx, idx)] = 1.0 / sigma;
        }
    }
    Ok(vt.transpose() * sigma_inv * u.transpose())
}

#[derive(Clone)]
struct SingularTermPlumbing {
    a: f64,
    s: DMatrix<f64>,
    d: DMatrix<f64>,
    b: DMatrix<f64>,
}

impl SingularTermPlumbing {
    fn new(singular_term: &faer_dense_mat, a: f64) -> Result<Self, String> {
        if singular_term.nrows() != singular_term.ncols() {
            return Err(format!(
                "singular term matrix must be square, got {}x{}",
                singular_term.nrows(),
                singular_term.ncols()
            ));
        }
        let s = faer_dense_to_dmatrix(singular_term);
        let pinv_s = pseudoinverse(&s)?;
        let identity = DMatrix::<f64>::identity(s.nrows(), s.ncols());
        let b = &identity - &(pinv_s * &s);
        let d = pseudoinverse(&(identity - &s))?;
        Ok(Self { a, s, d, b })
    }

    fn project_initial_guess(&self, y: &mut faer_dense_mat) {
        if y.ncols() == 0 {
            return;
        }
        let y_dense = faer_dense_to_dmatrix(y);
        let projected = &self.b * y_dense.column(0).clone_owned();
        for row in 0..y.nrows() {
            *y.get_mut(row, 0) = projected[row];
        }
    }

    fn wrap_rhs_output(
        &self,
        x: &faer_col,
        base: &faer_dense_mat,
        y: &faer_dense_mat,
    ) -> faer_dense_mat {
        let mut out = faer_dense_to_dmatrix(base);
        let y_dense = faer_dense_to_dmatrix(y);
        if x.nrows() == 0 {
            return dmatrix_to_faer_dense(&out);
        }

        if x[0] == self.a {
            let first = out.column(0).clone_owned();
            let transformed = &self.d * first;
            out.column_mut(0).copy_from(&transformed);

            for col in 1..x.nrows() {
                let denom = x[col] - self.a;
                if denom == 0.0 {
                    continue;
                }
                let correction = &self.s * y_dense.column(col).clone_owned();
                for row in 0..out.nrows() {
                    out[(row, col)] += correction[row] / denom;
                }
            }
        } else {
            for col in 0..x.nrows() {
                let denom = x[col] - self.a;
                if denom == 0.0 {
                    continue;
                }
                let correction = &self.s * y_dense.column(col).clone_owned();
                for row in 0..out.nrows() {
                    out[(row, col)] += correction[row] / denom;
                }
            }
        }

        dmatrix_to_faer_dense(&out)
    }
    #[allow(dead_code)]
    fn wrap_pointwise_jacobians(&self, x: &faer_col, jacobians: Vec<faer_mat>) -> Vec<faer_mat> {
        let mut wrapped = Vec::with_capacity(jacobians.len());
        for (idx, jacobian) in jacobians.into_iter().enumerate() {
            let mut dense = sparse_to_dmatrix(&jacobian);
            if x.nrows() > 0 && x[0] == self.a && idx == 0 {
                dense = &self.d * dense;
            } else {
                let denom = x[idx] - self.a;
                if denom != 0.0 {
                    dense += &self.s / denom;
                }
            }
            wrapped.push(dmatrix_to_sparse(&dense));
        }
        wrapped
    }
}

fn record_global_jacobian_diagnostics(
    calc_statistics: &mut HashMap<String, usize>,
    diagnostics: &GlobalJacobianDiagnostics,
) {
    calc_statistics.insert("global jacobian rows".to_string(), diagnostics.nrows);
    calc_statistics.insert("global jacobian cols".to_string(), diagnostics.ncols);
    calc_statistics.insert("global jacobian nnz".to_string(), diagnostics.nnz);
    calc_statistics.insert(
        "global jacobian dense equivalent bytes".to_string(),
        diagnostics.dense_equivalent_bytes,
    );
    calc_statistics.insert(
        "global jacobian sparse values bytes".to_string(),
        diagnostics.sparse_values_bytes,
    );
    calc_statistics.insert(
        "global jacobian sparse row index bytes".to_string(),
        diagnostics.sparse_row_index_bytes,
    );
    calc_statistics.insert(
        "global jacobian sparse col ptr bytes".to_string(),
        diagnostics.sparse_col_ptr_bytes,
    );
    calc_statistics.insert(
        "global jacobian sparse storage bytes".to_string(),
        diagnostics.sparse_storage_bytes,
    );
    calc_statistics.insert(
        "global jacobian dense equivalent kib".to_string(),
        diagnostics.dense_equivalent_kib,
    );
    calc_statistics.insert(
        "global jacobian dense equivalent mib".to_string(),
        diagnostics.dense_equivalent_mib,
    );
    calc_statistics.insert(
        "global jacobian sparse storage kib".to_string(),
        diagnostics.sparse_storage_kib,
    );
    calc_statistics.insert(
        "global jacobian sparse storage mib".to_string(),
        diagnostics.sparse_storage_mib,
    );
    calc_statistics.insert(
        "global jacobian dense to sparse permille".to_string(),
        diagnostics.dense_to_sparse_permille,
    );
    calc_statistics.insert(
        "global jacobian nonfinite values".to_string(),
        diagnostics.nonfinite_values,
    );
    calc_statistics.insert(
        "global jacobian structurally empty rows".to_string(),
        diagnostics.structurally_empty_rows,
    );
    calc_statistics.insert(
        "global jacobian structurally empty cols".to_string(),
        diagnostics.structurally_empty_cols,
    );
    calc_statistics.insert(
        "global jacobian numerically zero rows".to_string(),
        diagnostics.numerically_zero_rows,
    );
    calc_statistics.insert(
        "global jacobian numerically zero cols".to_string(),
        diagnostics.numerically_zero_cols,
    );
}
// demonstrating function or reminder function with main
#[allow(dead_code)]
pub fn test_faer_fn() {
    type faer_mat = SparseColMat<usize, f64>;
    type faer_col = Col<f64>; // Mat<f64>;
    // Example usage of faer types
    type faer_mat_sym = SymbolicSparseColMat<usize, f64>;
    // empty vector
    let mut vec1: faer_col = faer_col::zeros(2);
    // setting values of vector
    vec1[1] = 1.0;
    // 1st method to create sparse matrix
    // col_ptr: [0, 1, 1, 2, 2, 2]
    // - Column 0: 1 entry (row 0)
    // - Column 1: 0 entries
    // - Column 2: 1 entry (row 1)
    // - Column 3: 0 entries
    // - Column 4: 0 entries
    let nrows = 5;
    let ncols = 5;
    let col_ptr = vec![0usize, 1, 1, 2, 2, 2]; // length = ncols + 1
    let row_idx = vec![0usize, 1]; // length = col_ptr[ncols]
    let values = vec![1.0, 2.0]; // same length as row_idx
    // col_nnz in SymbolicSparseColMat (and related types) stands for column nonzero counts.
    // It is an optional vector that, if present, specifies for
    // each column the number of nonzero entries in that column.
    let symbolic = SymbolicSparseColMat::new_checked(nrows, ncols, col_ptr, None, row_idx);
    let matrix = SparseColMat::new(symbolic, values);
    // getting certain element from matrix
    let _element = matrix.get(1, 1).unwrap_or(&0.0);
    //creating sparse mattrix from vector of triplets
    let non_zero_triplet = vec![(0, 0, 1.0), (1, 1, 1.0)];
    let nrows = 2;
    let ncols = 2;

    let triplet: Vec<Triplet<usize, usize, f64>> = non_zero_triplet
        .iter()
        .map(|triplet| Triplet::new(triplet.0, triplet.1, triplet.2))
        .collect::<Vec<_>>();
    let new_matrix: SparseColMat<usize, f64> =
        SparseColMat::try_new_from_triplets(nrows, ncols, triplet.as_slice()).unwrap();
    // solving procedure for linear system
    let lhs: MatRef<f64> = vec1.as_mat();
    let LU = new_matrix.sp_lu().unwrap();
    let res: Mat<f64> = LU.solve(lhs);
    let res_vec: Vec<f64> = res.row_iter().map(|x| x[0]).collect();
    // calculation result
    let _res: faer_col = ColRef::from_slice(res_vec.as_slice()).to_owned();
}

/// Cache for the symbolic sparsity pattern of the global Jacobian matrix.
///
/// The global Jacobian has a fixed block structure determined by `n`, `m`, and `k`:
/// - `(m-1)` diagonal `n×n` blocks (identity + Jacobian contributions)
/// - `(m-1)` off-diagonal `n×n` blocks (identity + Jacobian contributions)
/// - BC row blocks: `(n+k)×n` for ya, `(n+k)×n` for yb
/// - Parameter column blocks (if `k > 0`): `(m-1)*n×k` and `(n+k)×k`
///
/// The cache stores the `SparseColMat` with the correct symbolic structure and
/// a mutable values buffer. On each Newton iteration, [`construct_global_jac_cached`]
/// writes numerical values directly into the buffer without rebuilding the sparsity
/// pattern or sorting triplets.
pub struct JacobianSparsityCache {
    matrix: SparseColMat<usize, f64>,
    n: usize,
    m: usize,
    k: usize,
}

impl JacobianSparsityCache {
    /// Build a new cache for the given dimensions.
    ///
    /// Constructs the full symbolic sparsity pattern (all possible non-zeros) and
    /// stores it as a `SparseColMat` with a zeroed values buffer.
    pub fn new(n: usize, m: usize, k: usize) -> Self {
        let matrix = Self::build_structure_triplets(n, m, k);
        JacobianSparsityCache { matrix, n, m, k }
    }

    /// Return a mutable reference to the values buffer.
    pub fn values_mut(&mut self) -> &mut [f64] {
        let (_, values) = self.matrix.parts_mut();
        values
    }

    /// Return an immutable reference to the cached matrix (for reading).
    pub fn matrix_ref(&self) -> &SparseColMat<usize, f64> {
        &self.matrix
    }

    /// Get the symbolic structure parts: (nrows, ncols, col_ptr, col_nnz, row_idx)
    pub fn symbolic_parts(&mut self) -> (usize, usize, &[usize], Option<&[usize]>, &[usize]) {
        let (sym_ref, _) = self.matrix.parts_mut();
        sym_ref.parts()
    }

    /// Return a clone of the underlying matrix (after values have been written).
    pub fn matrix_clone(&self) -> SparseColMat<usize, f64> {
        self.matrix.clone()
    }

    /// Build the full sparsity pattern as a `SparseColMat` by constructing triplets
    /// for every possible non-zero position in the global Jacobian.
    ///
    /// The pattern includes:
    /// - Identity entries on diagonal/off-diagonal block diagonals
    /// - Dense `n×n` blocks for `df_dy` / `df_dy_middle` contributions
    /// - BC blocks for `dbc_dya`, `dbc_dyb`, `dbc_dp`
    fn build_structure_triplets(n: usize, m: usize, k: usize) -> SparseColMat<usize, f64> {
        let total_size = (m - 1) * n + (n + k);
        let mut triplets = Vec::new();

        // Diagonal and off-diagonal blocks for each collocation interval
        for i in 0..(m - 1) {
            let row_start = i * n;

            // --- Diagonal block: column range [i*n, (i+1)*n) ---
            let col_start = i * n;

            // Identity part: -I on diagonal
            for idx in 0..n {
                triplets.push(Triplet::new(row_start + idx, col_start + idx, -1.0));
            }

            // Dense n×n block for -h/6 * df_dy[i]
            for r in 0..n {
                for c in 0..n {
                    triplets.push(Triplet::new(row_start + r, col_start + c, 0.0));
                }
            }

            // Dense n×n block for -h/3 * df_dy_middle[i]
            for r in 0..n {
                for c in 0..n {
                    triplets.push(Triplet::new(row_start + r, col_start + c, 0.0));
                }
            }

            // --- Off-diagonal block: column range [(i+1)*n, (i+2)*n) ---
            let col_start_off = (i + 1) * n;

            // Identity part: +I on diagonal
            for idx in 0..n {
                triplets.push(Triplet::new(row_start + idx, col_start_off + idx, 1.0));
            }

            // Dense n×n block for -h/6 * df_dy[i+1]
            for r in 0..n {
                for c in 0..n {
                    triplets.push(Triplet::new(row_start + r, col_start_off + c, 0.0));
                }
            }

            // Dense n×n block for -h/3 * df_dy_middle[i]
            for r in 0..n {
                for c in 0..n {
                    triplets.push(Triplet::new(row_start + r, col_start_off + c, 0.0));
                }
            }

            // Parameter column blocks (if k > 0)
            if k > 0 {
                let param_col_start = n * m;

                // -h/6 * df_dp[i]
                for r in 0..n {
                    for c in 0..k {
                        triplets.push(Triplet::new(row_start + r, param_col_start + c, 0.0));
                    }
                }

                // -h/6 * df_dp[i+1]
                for r in 0..n {
                    for c in 0..k {
                        triplets.push(Triplet::new(row_start + r, param_col_start + c, 0.0));
                    }
                }

                // -2h/3 * df_dp_middle[i]
                for r in 0..n {
                    for c in 0..k {
                        triplets.push(Triplet::new(row_start + r, param_col_start + c, 0.0));
                    }
                }
            }
        }

        // Boundary condition blocks
        let bc_row_start = (m - 1) * n;

        // dbc_dya block: (n+k) × n
        for r in 0..(n + k) {
            for c in 0..n {
                triplets.push(Triplet::new(bc_row_start + r, c, 0.0));
            }
        }

        // dbc_dyb block: (n+k) × n, starting at column (m-1)*n
        let yb_col_start = (m - 1) * n;
        for r in 0..(n + k) {
            for c in 0..n {
                triplets.push(Triplet::new(bc_row_start + r, yb_col_start + c, 0.0));
            }
        }

        // dbc_dp block: (n+k) × k (if k > 0)
        if k > 0 {
            let param_col_start = n * m;
            for r in 0..(n + k) {
                for c in 0..k {
                    triplets.push(Triplet::new(bc_row_start + r, param_col_start + c, 0.0));
                }
            }
        }

        SparseColMat::try_new_from_triplets(total_size, total_size, &triplets)
            .expect("Failed to create sparsity cache matrix")
    }
}

/// Machine epsilon for floating point arithmetic
const EPS: f64 = f64::EPSILON;

/// Result structure for BVP solver
#[derive(Clone)]
pub struct BVPResult {
    /// Solution as cubic spline interpolator (PPoly)
    pub sol: Option<PPoly>,
    /// Found parameters (if any)
    pub p: Option<faer_col>,
    /// Final mesh nodes
    pub x: faer_col,
    /// Solution values at mesh nodes
    pub y: faer_dense_mat,
    /// Solution derivatives at mesh nodes  
    pub yp: faer_dense_mat,
    /// RMS residuals over each mesh interval
    pub rms_residuals: faer_col,
    /// Number of iterations
    pub niter: usize,
    /// Status code (0=success, 1=max_nodes, 2=singular, 3=bc_tol)
    pub status: i32,
    /// Status message
    pub message: String,
    /// Success flag
    pub success: bool,
    //
    pub calc_statistics: HashMap<String, usize>,
    //
    pub custom_timer: CustomTimer,
    //
    pub strategy_params: Option<HashMap<String, Vec<f64>>>,
}
impl Default for BVPResult {
    fn default() -> Self {
        // Initialize with default values
        let vec_of_tuples = vec![
            ("number of iterations".to_string(), 0),
            ("number of solving linear systems".to_string(), 0),
            ("number of jacobians recalculations".to_string(), 0),
            ("number of grid refinements".to_string(), 0),
        ];
        let Hashmap_statistics: HashMap<String, usize> = vec_of_tuples.into_iter().collect();
        let max_iter = 8;
        let max_njev = 4;
        let sigma = 0.2; // Armijo constant
        let tau = 0.5; // Step size decrease factor
        let n_trial = 4; // Max backtracking step
        let strategy_params: Option<HashMap<String, Vec<f64>>> = Some(
            vec![
                ("max_iter".to_string(), vec![max_iter as f64]),
                ("max_njev".to_string(), vec![max_njev as f64]),
                ("sigma".to_string(), vec![sigma]),
                ("tau".to_string(), vec![tau]),
                ("n_trial".to_string(), vec![n_trial as f64]),
            ]
            .into_iter()
            .collect(),
        );
        BVPResult {
            sol: None,
            p: None,
            x: faer_col::zeros(0),
            y: faer_dense_mat::zeros(0, 0),
            yp: faer_dense_mat::zeros(0, 0),
            rms_residuals: faer_col::zeros(0),
            niter: 0,
            status: 0,
            message: "".to_string(),
            success: false,
            calc_statistics: Hashmap_statistics,
            custom_timer: CustomTimer::new(),
            strategy_params: strategy_params,
        }
    }
}

impl BVPResult {
    /// Set strategy parameters for the BVP solver
    pub fn set_strategy_params(&mut self, params: HashMap<String, Vec<f64>>) {
        self.strategy_params = Some(params);
    }
}

fn default_strategy_params() -> HashMap<String, Vec<f64>> {
    BVPResult::default()
        .strategy_params
        .expect("BVPResult default strategy parameters should exist")
}

fn strategy_param_f64(
    strategy_params: Option<&HashMap<String, Vec<f64>>>,
    key: &str,
    fallback: f64,
) -> f64 {
    strategy_params
        .and_then(|params| params.get(key))
        .and_then(|values| values.first().copied())
        .unwrap_or(fallback)
}

fn strategy_param_usize(
    strategy_params: Option<&HashMap<String, Vec<f64>>>,
    key: &str,
    fallback: usize,
) -> usize {
    strategy_param_f64(strategy_params, key, fallback as f64) as usize
}
/// Function type for ODE right-hand side evaluation
/// Arguments: (x, y, p) where x is scalar, y is n-dimensional, p is k-dimensional parameters
pub type ODEFunction = dyn Fn(&faer_col, &faer_dense_mat, &faer_col) -> faer_dense_mat;

/// Function type for boundary condition evaluation  
/// Arguments: (ya, yb, p) where ya, yb are n-dimensional boundary values, p is parameters
/// Returns: (n+k)-dimensional boundary condition residuals
pub type BCFunction = dyn Fn(&faer_col, &faer_col, &faer_col) -> faer_col;

/// Function type for ODE Jacobian evaluation (optional)
/// Returns: (df_dy, df_dp) where df_dy is (n,n,m), df_dp is (n,k,m) or None
pub type ODEJacobian =
    dyn Fn(&faer_col, &faer_dense_mat, &faer_col) -> (Vec<faer_mat>, Option<Vec<faer_mat>>);

/// Function type for boundary condition Jacobian evaluation (optional)  
/// Returns: (dbc_dya, dbc_dyb, dbc_dp) where each is appropriately sized matrix or None
pub type BCJacobian =
    dyn Fn(&faer_col, &faer_col, &faer_col) -> (faer_mat, faer_mat, Option<faer_mat>);

/// Estimate derivatives of an ODE system RHS with forward differences
///
/// # Arguments
/// * `fun` - ODE function to differentiate
/// * `x` - Mesh points (m,)
/// * `y` - Solution values at mesh points (n, m)  
/// * `p` - Parameters (k,)
/// * `f0` - Pre-computed function values (optional)
///
/// # Returns
/// * `df_dy` - Derivatives w.r.t. y: Vec of (n,n) matrices, one per mesh point
/// * `df_dp` - Derivatives w.r.t. p: Vec of (n,k) matrices, one per mesh point (or None)
pub fn estimate_fun_jac(
    fun: &ODEFunction,
    x: &faer_col,
    y: &faer_dense_mat,
    p: &faer_col,
    f0: Option<&faer_dense_mat>,
) -> (Vec<faer_mat>, Option<Vec<faer_mat>>) {
    let (n, m) = (y.nrows(), y.ncols());

    let f0_computed;
    let f0_ref = match f0 {
        Some(f) => f,
        None => {
            f0_computed = fun(x, y, p);
            &f0_computed
        }
    };

    let mut df_dy = Vec::with_capacity(m);

    // Compute df/dy for each mesh point
    for col in 0..m {
        let mut triplets = Vec::new();

        for i in 0..n {
            let mut y_perturbed = y.clone();
            let h = EPS.sqrt() * (1.0 + y.get(i, col).abs());
            *y_perturbed.get_mut(i, col) += h;

            let x_slice = x.clone();
            let f_new = fun(&x_slice, &y_perturbed, p);

            for row in 0..n {
                let val = (f_new.get(row, col) - f0_ref.get(row, col)) / h;
                if val.abs() > 1e-15 {
                    triplets.push(Triplet::new(row, i, val));
                }
            }
        }
        let jacobian = SparseColMat::try_new_from_triplets(n, n, &triplets).unwrap();
        df_dy.push(jacobian);
    }

    // Compute df/dp if parameters exist
    let df_dp = if p.nrows() == 0 {
        None
    } else {
        let k = p.nrows();
        let mut df_dp_vec = Vec::with_capacity(m);

        for col in 0..m {
            let mut triplets = Vec::new();

            for i in 0..k {
                let mut p_perturbed = p.clone();
                let h = EPS.sqrt() * (1.0 + p[i].abs());
                p_perturbed[i] += h;

                let f_new = fun(x, y, &p_perturbed);

                for row in 0..n {
                    let val = (f_new.get(row, col) - f0_ref.get(row, col)) / h;
                    if val.abs() > 1e-15 {
                        triplets.push(Triplet::new(row, i, val));
                    }
                }
            }
            let param_jacobian = SparseColMat::try_new_from_triplets(n, k, &triplets).unwrap();
            df_dp_vec.push(param_jacobian);
        }
        Some(df_dp_vec)
    };

    (df_dy, df_dp)
}

/// Estimate derivatives of boundary conditions with forward differences
///
/// # Arguments  
/// * `bc` - Boundary condition function
/// * `ya` - Left boundary values (n,)
/// * `yb` - Right boundary values (n,)
/// * `p` - Parameters (k,)
/// * `bc0` - Pre-computed boundary condition values (optional)
///
/// # Returns
/// * `dbc_dya` - Derivatives w.r.t. ya: (n+k, n) matrix
/// * `dbc_dyb` - Derivatives w.r.t. yb: (n+k, n) matrix  
/// * `dbc_dp` - Derivatives w.r.t. p: (n+k, k) matrix (or None)
pub fn estimate_bc_jac(
    bc: &BCFunction,
    ya: &faer_col,
    yb: &faer_col,
    p: &faer_col,
    bc0: Option<&faer_col>,
) -> (faer_mat, faer_mat, Option<faer_mat>) {
    let n = ya.nrows();
    let k = p.nrows();

    let bc0_computed;
    let bc0_ref = match bc0 {
        Some(bc_val) => bc_val,
        None => {
            bc0_computed = bc(ya, yb, p);
            &bc0_computed
        }
    };

    // Compute dbc/dya
    let mut triplets_ya = Vec::new();
    for i in 0..n {
        let mut ya_perturbed = ya.clone();
        let h = EPS.sqrt() * (1.0 + ya[i].abs());
        ya_perturbed[i] += h;

        let bc_new = bc(&ya_perturbed, yb, p);
        for row in 0..(n + k) {
            let val = (bc_new[row] - bc0_ref[row]) / h;
            if val.abs() > 1e-15 {
                triplets_ya.push(Triplet::new(row, i, val));
            }
        }
    }
    let dbc_dya = SparseColMat::try_new_from_triplets(n + k, n, &triplets_ya).unwrap();

    // Compute dbc/dyb
    let mut triplets_yb = Vec::new();
    for i in 0..n {
        let mut yb_perturbed = yb.clone();
        let h = EPS.sqrt() * (1.0 + yb[i].abs());
        yb_perturbed[i] += h;

        let bc_new = bc(ya, &yb_perturbed, p);
        for row in 0..(n + k) {
            let val = (bc_new[row] - bc0_ref[row]) / h;
            if val.abs() > 1e-15 {
                triplets_yb.push(Triplet::new(row, i, val));
            }
        }
    }
    let dbc_dyb = SparseColMat::try_new_from_triplets(n + k, n, &triplets_yb).unwrap();

    // Compute dbc/dp if parameters exist
    let dbc_dp = if k == 0 {
        None
    } else {
        let mut triplets_p = Vec::new();
        for i in 0..k {
            let mut p_perturbed = p.clone();
            let h = EPS.sqrt() * (1.0 + p[i].abs());
            p_perturbed[i] += h;

            let bc_new = bc(ya, yb, &p_perturbed);
            for row in 0..(n + k) {
                let val = (bc_new[row] - bc0_ref[row]) / h;
                if val.abs() > 1e-15 {
                    triplets_p.push(Triplet::new(row, i, val));
                }
            }
        }
        Some(SparseColMat::try_new_from_triplets(n + k, k, &triplets_p).unwrap())
    };

    (dbc_dya, dbc_dyb, dbc_dp)
}

/// Compute indices for the collocation system Jacobian construction
///
/// # Arguments
/// * `n` - Number of equations in ODE system
/// * `m` - Number of mesh nodes  
/// * `k` - Number of unknown parameters
///
/// # Returns
/// * `i_indices` - Row indices for sparse matrix construction
/// * `j_indices` - Column indices for sparse matrix construction  
pub fn compute_jac_indices(n: usize, m: usize, k: usize) -> (Vec<usize>, Vec<usize>) {
    // Pattern matches the Python implementation structure
    let mut i_indices = Vec::new();
    let mut j_indices = Vec::new();

    // Block 1: m-1 diagonal n×n blocks for collocation residuals
    for block in 0..(m - 1) {
        for i in 0..n {
            for j in 0..n {
                i_indices.push(block * n + i);
                j_indices.push(block * n + j);
            }
        }
    }

    // Block 2: m-1 off-diagonal n×n blocks for collocation residuals
    for block in 0..(m - 1) {
        for i in 0..n {
            for j in 0..n {
                i_indices.push(block * n + i);
                j_indices.push((block + 1) * n + j);
            }
        }
    }

    // Block 3: (n+k)×n block for dependency of BC on ya
    for i in 0..(n + k) {
        for j in 0..n {
            i_indices.push((m - 1) * n + i);
            j_indices.push(j);
        }
    }

    // Block 4: (n+k)×n block for dependency of BC on yb
    for i in 0..(n + k) {
        for j in 0..n {
            i_indices.push((m - 1) * n + i);
            j_indices.push((m - 1) * n + j);
        }
    }

    // Block 5: (m-1)*n×k block for dependency of collocation on p
    if k > 0 {
        for block in 0..(m - 1) {
            for i in 0..n {
                for j in 0..k {
                    i_indices.push(block * n + i);
                    j_indices.push(m * n + j);
                }
            }
        }

        // Block 6: (n+k)×k block for dependency of BC on p
        for i in 0..(n + k) {
            for j in 0..k {
                i_indices.push((m - 1) * n + i);
                j_indices.push(m * n + j);
            }
        }
    }

    (i_indices, j_indices)
}

/// Evaluate collocation residuals
///
/// This function implements the core collocation method. The solution is sought
/// as a cubic C1 continuous spline with derivatives matching the ODE RHS at given nodes.
/// Collocation conditions are formed from equality of spline derivatives and RHS
/// of the ODE system at middle points between nodes.
///
/// # Arguments
/// * `fun` - ODE function f(x, y, p)
/// * `y` - Solution values at mesh nodes (n, m)
/// * `p` - Parameters (k,)  
/// * `x` - Mesh nodes (m,)
/// * `h` - Mesh intervals (m-1,)
///
/// # Returns
/// * `col_res` - Collocation residuals at middle points (n, m-1)
/// * `y_middle` - Spline values at middle points (n, m-1)
/// * `f` - RHS values at mesh nodes (n, m)
/// * `f_middle` - RHS values at middle points (n, m-1)
pub fn collocation_fun(
    fun: &ODEFunction,
    y: &faer_dense_mat,
    p: &faer_col,
    x: &faer_col,
    h: &faer_col,
) -> (
    faer_dense_mat,
    faer_dense_mat,
    faer_dense_mat,
    faer_dense_mat,
) {
    let (n, m) = (y.nrows(), y.ncols());

    // Evaluate RHS at mesh nodes
    let f = fun(x, y, p);

    // Compute solution values at middle points using cubic interpolation formula
    let mut y_middle = faer_dense_mat::zeros(n, m - 1);
    for i in 0..n {
        for j in 0..(m - 1) {
            *y_middle.get_mut(i, j) = 0.5 * (y.get(i, j + 1) + y.get(i, j))
                - 0.125 * h[j] * (f.get(i, j + 1) - f.get(i, j));
        }
    }

    // Evaluate RHS at middle points
    // The ODE function already supports batched evaluation across multiple points
    // (same as the mesh-node call above). Build a single column vector of all
    // midpoints and evaluate in one call for optimal performance.
    let mut x_mid_all = faer_col::zeros(m - 1);
    for j in 0..(m - 1) {
        x_mid_all[j] = x[j] + 0.5 * h[j];
    }
    let f_middle = fun(&x_mid_all, &y_middle, p);

    // Compute collocation residuals
    let mut col_res = faer_dense_mat::zeros(n, m - 1);
    for i in 0..n {
        for j in 0..(m - 1) {
            *col_res.get_mut(i, j) = y.get(i, j + 1)
                - y.get(i, j)
                - h[j] / 6.0 * (f.get(i, j) + f.get(i, j + 1) + 4.0 * f_middle.get(i, j));
        }
    }

    (col_res, y_middle, f, f_middle)
}

/// Construct the Jacobian of the collocation system using sparse matrices
pub fn construct_global_jac(
    n: usize,
    m: usize,
    k: usize,
    h: &faer_col,
    df_dy: &[faer_mat],
    df_dy_middle: &[faer_mat],
    df_dp: Option<&[faer_mat]>,
    df_dp_middle: Option<&[faer_mat]>,
    dbc_dya: &faer_mat,
    dbc_dyb: &faer_mat,
    dbc_dp: Option<&faer_mat>,
) -> faer_mat {
    let total_size = (m - 1) * n + (n + k);
    let mut triplets = Vec::new();

    // Process diagonal and off-diagonal blocks for collocation residuals
    for i in 0..(m - 1) {
        let h_i = h[i];

        // For sparse matrices, we'll compute the blocks directly as triplets
        // Diagonal block: -I - h/6*(df_dy[i] + 2*df_dy_middle[i]) - h²/12*df_dy_middle[i]*df_dy[i]

        // Add diagonal block entries
        let row_start = i * n;
        let col_start = i * n;

        // Identity part: -I
        for idx in 0..n {
            triplets.push(Triplet::new(row_start + idx, col_start + idx, -1.0));
        }

        // Add -h/6 * df_dy[i]
        for triplt in df_dy[i].triplet_iter() {
            let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
            let coeff = -h_i / 6.0 * val;
            if coeff.abs() > 1e-15 {
                triplets.push(Triplet::new(
                    row_start + row_idx,
                    col_start + col_idx,
                    coeff,
                ));
            }
        }

        // Add -h/3 * df_dy_middle[i] (2 * h/6)
        for triplt in df_dy_middle[i].triplet_iter() {
            let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
            let coeff = -h_i / 3.0 * val;
            if coeff.abs() > 1e-15 {
                triplets.push(Triplet::new(
                    row_start + row_idx,
                    col_start + col_idx,
                    coeff,
                ));
            }
        }

        // Off-diagonal block: I - h/6*(df_dy[i+1] + 2*df_dy_middle[i]) + h²/12*df_dy_middle[i]*df_dy[i+1]

        // Add off-diagonal block entries
        let col_start_off = (i + 1) * n;

        // Identity part: I
        for idx in 0..n {
            triplets.push(Triplet::new(row_start + idx, col_start_off + idx, 1.0));
        }

        // Add -h/6 * df_dy[i+1]
        for triplt in df_dy[i + 1].triplet_iter() {
            let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
            let coeff = -h_i / 6.0 * val;
            if coeff.abs() > 1e-15 {
                triplets.push(Triplet::new(
                    row_start + row_idx,
                    col_start_off + col_idx,
                    coeff,
                ));
            }
        }

        // Add -h/3 * df_dy_middle[i]
        for triplt in df_dy_middle[i].triplet_iter() {
            let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
            let coeff = -h_i / 3.0 * val;
            if coeff.abs() > 1e-15 {
                triplets.push(Triplet::new(
                    row_start + row_idx,
                    col_start_off + col_idx,
                    coeff,
                ));
            }
        }

        // Handle parameter derivatives if present
        if let (Some(df_dp_vec), Some(df_dp_middle_vec)) = (df_dp, df_dp_middle) {
            let param_col_start = n * m;

            // Add -h/6 * df_dp[i]
            for triplt in df_dp_vec[i].triplet_iter() {
                let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
                let coeff = -h_i / 6.0 * val;
                if coeff.abs() > 1e-15 {
                    triplets.push(Triplet::new(
                        row_start + row_idx,
                        param_col_start + col_idx,
                        coeff,
                    ));
                }
            }

            // Add -h/6 * df_dp[i+1]
            for triplt in df_dp_vec[i + 1].triplet_iter() {
                let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
                let coeff = -h_i / 6.0 * val;
                if coeff.abs() > 1e-15 {
                    triplets.push(Triplet::new(
                        row_start + row_idx,
                        param_col_start + col_idx,
                        coeff,
                    ));
                }
            }

            // Add -2h/3 * df_dp_middle[i]
            for triplt in df_dp_middle_vec[i].triplet_iter() {
                let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
                let coeff = -2.0 * h_i / 3.0 * val;
                if coeff.abs() > 1e-15 {
                    triplets.push(Triplet::new(
                        row_start + row_idx,
                        param_col_start + col_idx,
                        coeff,
                    ));
                }
            }
        }
    }

    // Insert boundary condition blocks
    let bc_row_start = (m - 1) * n;

    // dbc_dya block - dependency on ya
    for triplt in dbc_dya.triplet_iter() {
        let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
        if val.abs() > 1e-15 {
            triplets.push(Triplet::new(bc_row_start + row_idx, col_idx, *val));
        }
    }

    // dbc_dyb block - dependency on yb
    let yb_col_start = (m - 1) * n;
    for triplt in dbc_dyb.triplet_iter() {
        let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
        if val.abs() > 1e-15 {
            triplets.push(Triplet::new(
                bc_row_start + row_idx,
                yb_col_start + col_idx,
                *val,
            ));
        }
    }

    // dbc_dp block - dependency on parameters
    if let Some(dbc_dp_mat) = dbc_dp {
        let param_col_start = n * m;
        for triplt in dbc_dp_mat.triplet_iter() {
            let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
            if val.abs() > 1e-15 {
                triplets.push(Triplet::new(
                    bc_row_start + row_idx,
                    param_col_start + col_idx,
                    *val,
                ));
            }
        }
    }

    // Create sparse matrix from triplets
    SparseColMat::try_new_from_triplets(total_size, total_size, &triplets)
        .expect("Failed to create sparse matrix")
}

/// Construct the global Jacobian matrix using a pre-computed sparsity cache.
///
/// This is the cached counterpart of [`construct_global_jac`]. Instead of building
/// triplets and sorting them (which is O(nnz log nnz)), it iterates over the cache's
/// CSC structure directly and writes numerical values into the pre-allocated buffer.
///
/// The cache's symbolic pattern includes ALL possible non-zeros (identity entries +
/// dense n×n blocks for df_dy/df_dy_middle + BC blocks). For each (col, row) pair
/// in the cache's column pointers and row indices, this function determines which
/// block the entry belongs to and computes the numerical value by **accumulating**
/// all contributions (identity + Jacobian) at that position.
///
/// # Correctness: Accumulation, not replacement
///
/// This function uses a **faithful** approach: it builds the same triplets as
/// [`construct_global_jac`], sorts them by (col, row) (matching `try_new_from_triplets`
/// internal behavior), then walks through the cache's CSC structure and sorted triplets
/// simultaneously, accumulating values at matching positions. This guarantees bit-identical
/// results with the original function.
pub fn construct_global_jac_cached(
    n: usize,
    m: usize,
    k: usize,
    h: &faer_col,
    df_dy: &[faer_mat],
    df_dy_middle: &[faer_mat],
    df_dp: Option<&[faer_mat]>,
    df_dp_middle: Option<&[faer_mat]>,
    dbc_dya: &faer_mat,
    dbc_dyb: &faer_mat,
    dbc_dp: Option<&faer_mat>,
    cache: &mut JacobianSparsityCache,
) {
    // Step 1: Build triplets using EXACTLY the same logic as construct_global_jac
    let _total_size = (m - 1) * n + (n + k);
    let mut triplets = Vec::new();

    // Process diagonal and off-diagonal blocks for collocation residuals
    for i in 0..(m - 1) {
        let h_i = h[i];

        // Diagonal block: -I - h/6*(df_dy[i] + 2*df_dy_middle[i])
        let row_start = i * n;
        let col_start = i * n;

        // Identity part: -I
        for idx in 0..n {
            triplets.push(Triplet::new(row_start + idx, col_start + idx, -1.0));
        }

        // Add -h/6 * df_dy[i]
        for triplt in df_dy[i].triplet_iter() {
            let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
            let coeff = -h_i / 6.0 * val;
            if coeff.abs() > 1e-15 {
                triplets.push(Triplet::new(
                    row_start + row_idx,
                    col_start + col_idx,
                    coeff,
                ));
            }
        }

        // Add -h/3 * df_dy_middle[i]
        for triplt in df_dy_middle[i].triplet_iter() {
            let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
            let coeff = -h_i / 3.0 * val;
            if coeff.abs() > 1e-15 {
                triplets.push(Triplet::new(
                    row_start + row_idx,
                    col_start + col_idx,
                    coeff,
                ));
            }
        }

        // Off-diagonal block: I - h/6*(df_dy[i+1] + 2*df_dy_middle[i])
        let col_start_off = (i + 1) * n;

        // Identity part: I
        for idx in 0..n {
            triplets.push(Triplet::new(row_start + idx, col_start_off + idx, 1.0));
        }

        // Add -h/6 * df_dy[i+1]
        for triplt in df_dy[i + 1].triplet_iter() {
            let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
            let coeff = -h_i / 6.0 * val;
            if coeff.abs() > 1e-15 {
                triplets.push(Triplet::new(
                    row_start + row_idx,
                    col_start_off + col_idx,
                    coeff,
                ));
            }
        }

        // Add -h/3 * df_dy_middle[i]
        for triplt in df_dy_middle[i].triplet_iter() {
            let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
            let coeff = -h_i / 3.0 * val;
            if coeff.abs() > 1e-15 {
                triplets.push(Triplet::new(
                    row_start + row_idx,
                    col_start_off + col_idx,
                    coeff,
                ));
            }
        }

        // Handle parameter derivatives if present
        if let (Some(df_dp_vec), Some(df_dp_middle_vec)) = (df_dp, df_dp_middle) {
            let param_col_start = n * m;

            // Add -h/6 * df_dp[i]
            for triplt in df_dp_vec[i].triplet_iter() {
                let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
                let coeff = -h_i / 6.0 * val;
                if coeff.abs() > 1e-15 {
                    triplets.push(Triplet::new(
                        row_start + row_idx,
                        param_col_start + col_idx,
                        coeff,
                    ));
                }
            }

            // Add -h/6 * df_dp[i+1]
            for triplt in df_dp_vec[i + 1].triplet_iter() {
                let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
                let coeff = -h_i / 6.0 * val;
                if coeff.abs() > 1e-15 {
                    triplets.push(Triplet::new(
                        row_start + row_idx,
                        param_col_start + col_idx,
                        coeff,
                    ));
                }
            }

            // Add -2h/3 * df_dp_middle[i]
            for triplt in df_dp_middle_vec[i].triplet_iter() {
                let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
                let coeff = -2.0 * h_i / 3.0 * val;
                if coeff.abs() > 1e-15 {
                    triplets.push(Triplet::new(
                        row_start + row_idx,
                        param_col_start + col_idx,
                        coeff,
                    ));
                }
            }
        }
    }

    // Insert boundary condition blocks
    let bc_row_start = (m - 1) * n;

    // dbc_dya block - dependency on ya
    for triplt in dbc_dya.triplet_iter() {
        let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
        if val.abs() > 1e-15 {
            triplets.push(Triplet::new(bc_row_start + row_idx, col_idx, *val));
        }
    }

    // dbc_dyb block - dependency on yb
    let yb_col_start = (m - 1) * n;
    for triplt in dbc_dyb.triplet_iter() {
        let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
        if val.abs() > 1e-15 {
            triplets.push(Triplet::new(
                bc_row_start + row_idx,
                yb_col_start + col_idx,
                *val,
            ));
        }
    }

    // dbc_dp block - dependency on parameters
    if let Some(dbc_dp_mat) = dbc_dp {
        let param_col_start = n * m;
        for triplt in dbc_dp_mat.triplet_iter() {
            let (row_idx, col_idx, val) = (triplt.row, triplt.col, triplt.val);
            if val.abs() > 1e-15 {
                triplets.push(Triplet::new(
                    bc_row_start + row_idx,
                    param_col_start + col_idx,
                    *val,
                ));
            }
        }
    }

    // Step 2: Sort triplets by (col, row) — matching try_new_from_triplets internal behavior
    triplets.sort_by(|a, b| {
        if a.col != b.col {
            a.col.cmp(&b.col)
        } else {
            a.row.cmp(&b.row)
        }
    });

    // Step 3: Zero out the cache's values buffer
    let (sym_ref, values) = cache.matrix.parts_mut();
    let (_nrows, _ncols, col_ptr, col_nnz, row_idx) = sym_ref.parts();
    for v in values.iter_mut() {
        *v = 0.0;
    }

    // Step 4: Walk through cache's CSC structure and sorted triplets simultaneously,
    //         accumulating values at matching positions.
    // The cache's CSC structure has entries for ALL possible non-zeros (including those
    // that would be filtered by coeff.abs() > 1e-15). The triplets only contain entries
    // that passed the filter. So we iterate over the cache's CSC positions and look for
    // matching triplets.
    let mut triplet_idx = 0;
    for col in 0.._ncols {
        let start = col_ptr[col];
        let nnz = match col_nnz {
            Some(cnts) => cnts[col],
            None => col_ptr[col + 1] - col_ptr[col],
        };
        let end = start + nnz;

        for pos in start..end {
            let row = row_idx[pos];

            // Accumulate all triplets at this (col, row) position
            while triplet_idx < triplets.len()
                && (triplets[triplet_idx].col < col
                    || (triplets[triplet_idx].col == col && triplets[triplet_idx].row < row))
            {
                triplet_idx += 1;
            }

            // Sum all triplets at this exact (col, row)
            while triplet_idx < triplets.len()
                && triplets[triplet_idx].col == col
                && triplets[triplet_idx].row == row
            {
                values[pos] += triplets[triplet_idx].val;
                triplet_idx += 1;
            }
        }
    }
}

fn record_linear_backend_route(
    calc_statistics: &mut HashMap<String, usize>,
    profile: &crate::numerical::BVP_sci::BVP_sci_bordered_banded::BvpSciBorderedBandedProfile,
) {
    let route_key = match profile.recommended_route {
        BvpSciBandedRoute::FullScalarBanded => "bvp sci route full scalar banded",
        BvpSciBandedRoute::BorderedBanded => "bvp sci route bordered banded candidate",
        BvpSciBandedRoute::SparseFallback => "bvp sci route sparse fallback",
    };
    *calc_statistics.entry(route_key.to_string()).or_insert(0) += 1;
    calc_statistics.insert("bvp sci full scalar kl".to_string(), profile.full_scalar.kl);
    calc_statistics.insert("bvp sci full scalar ku".to_string(), profile.full_scalar.ku);
    calc_statistics.insert(
        "bvp sci collocation scalar kl".to_string(),
        profile.collocation_scalar.kl,
    );
    calc_statistics.insert(
        "bvp sci collocation scalar ku".to_string(),
        profile.collocation_scalar.ku,
    );
    if let Some(amp) = profile.full_amplification {
        calc_statistics.insert(
            "bvp sci full scalar amp milli".to_string(),
            (amp * 1000.0).round() as usize,
        );
    }
    if let Some(amp) = profile.collocation_amplification {
        calc_statistics.insert(
            "bvp sci collocation scalar amp milli".to_string(),
            (amp * 1000.0).round() as usize,
        );
    }
}

fn add_elapsed_micros_stat(
    calc_statistics: &mut HashMap<String, usize>,
    key: &str,
    started_at: Instant,
) {
    let micros = started_at.elapsed().as_micros();
    let micros = usize::try_from(micros).unwrap_or(usize::MAX);
    let entry = calc_statistics.entry(key.to_string()).or_insert(0);
    *entry = entry.saturating_add(micros);
}

fn factor_banded_and_solve(
    jac_matrix: &faer_mat,
    rhs: &faer_col,
) -> Result<(LapackStyleBandedLuFaithful, faer_col), crate::somelinalg::banded::BandedError> {
    let banded = sparse_global_jac_to_banded(jac_matrix)?;
    let mut lu = LapackStyleBandedLuFaithful::new(banded.n(), banded.kl(), banded.ku())?;
    lu.factor_from(&banded)?;
    let step = solve_with_banded_factor(&lu, rhs)?;
    Ok((lu, step))
}

fn solve_with_banded_factor(
    lu: &LapackStyleBandedLuFaithful,
    rhs: &faer_col,
) -> Result<faer_col, crate::somelinalg::banded::BandedError> {
    let mut values = (0..rhs.nrows()).map(|i| rhs[i]).collect::<Vec<_>>();
    lu.solve_in_place(&mut values)?;
    Ok(faer_col::from_fn(values.len(), |i| values[i]))
}

/// Solve the nonlinear collocation system by Newton's method with sparse matrices
///
/// This version uses faer's sparse matrix capabilities for better performance.
///
/// # Arguments
/// * `n` - Number of equations in the ODE system
/// * `m` - Number of nodes in the mesh
/// * `h` - Mesh intervals (m-1,)
/// * `fun` - ODE function
/// * `bc` - Boundary condition function
/// * `fun_jac` - ODE Jacobian function (optional)
/// * `bc_jac` - BC Jacobian function (optional)
/// * `y` - Initial guess for solution values at mesh nodes (n, m)
/// * `p` - Initial guess for parameters (k,)
/// * `x` - Mesh nodes (m,)
/// * `bvp_tol` - BVP tolerance
/// * `bc_tol` - Boundary condition tolerance
///
/// # Returns
/// * `y` - Final solution values at mesh nodes
/// * `p` - Final parameter values  
/// * `singular` - True if Jacobian was singular
pub fn solve_newton(
    n: usize,
    m: usize,
    h: &faer_col,
    fun: &ODEFunction,
    bc: &BCFunction,
    fun_jac: Option<&ODEJacobian>,
    banded_fun_jac: Option<&ODEBandedJacobian>,
    bc_jac: Option<&BCJacobian>,
    mut y: faer_dense_mat,
    mut p: faer_col,
    x: &faer_col,
    bvp_tol: f64,
    bc_tol: f64,
    custom_timer: &mut CustomTimer,
    calc_statistics: &mut HashMap<String, usize>,
    mut sparsity_cache: Option<&mut JacobianSparsityCache>,
    strategy_params: Option<&HashMap<String, Vec<f64>>>,
    linear_solve_policy: BvpSciLinearSolvePolicy,
) -> (faer_dense_mat, faer_col, bool) {
    let k = p.nrows();

    // Newton iteration parameters
    let defaults = default_strategy_params();
    let strategy_params = strategy_params.or(Some(&defaults));
    let max_iter = strategy_param_usize(strategy_params, "max_iter", 8);
    let max_njev = strategy_param_usize(strategy_params, "max_njev", 4);
    let sigma = strategy_param_f64(strategy_params, "sigma", 0.2);
    let tau = strategy_param_f64(strategy_params, "tau", 0.5);
    let n_trial = strategy_param_usize(strategy_params, "n_trial", 4);

    // Tolerance for collocation residuals
    let tol_r: faer_col = faer_col::from_fn(h.nrows(), |i| 2.0 / 3.0 * h[i] * 5e-2 * bvp_tol);

    let mut njev = 0;
    let mut singular = false;
    let mut recompute_jac = true;
    let mut lu_solver: Option<faer::sparse::linalg::solvers::Lu<_, _>> = None;
    let mut banded_solver: Option<LapackStyleBandedLuFaithful> = None;
    let mut bordered_factorization: Option<BvpSciBorderedStructuredFactorization> = None;
    let mut bordered_workspace: Option<BvpSciBorderedSolveWorkspace> = None;
    let mut cost = 0.0;

    info!("\n \n Starting Newton iterations for collocation system \n \n");

    for iteration in 0..max_iter {
        info!("\n \n Newton iteration {} \n \n", iteration + 1);
        custom_timer.fun_tic();
        // Compute collocation residuals and function values
        let (col_res, y_middle, f, f_middle) = collocation_fun(fun, &y, &p, x, h);
        custom_timer.fun_tac();
        // Extract first and last columns for boundary conditions
        let ya = faer_col::from_fn(n, |i| *y.get(i, 0));
        let yb = faer_col::from_fn(n, |i| *y.get(i, m - 1));
        let bc_res = bc(&ya, &yb, &p);

        // Stack residuals into single vector (Python uses Fortran order: column-major)
        let mut res: Col<f64> = faer_col::zeros((m - 1) * n + (n + k));

        // Collocation residuals (column-major order like Python)
        for j in 0..(m - 1) {
            for i in 0..n {
                res[j * n + i] = *col_res.get(i, j);
            }
        }

        // Boundary condition residuals
        for i in 0..(n + k) {
            res[(m - 1) * n + i] = bc_res[i];
        }

        if recompute_jac {
            info!("\n \n JACOBIAN (RE)CALCULATED! \n \n");
            custom_timer.jac_tic();
            *calc_statistics
                .entry("number of jacobians recalculations".to_string())
                .or_insert(0) += 1;
            if let Some(banded_jac_fn) = banded_fun_jac {
                if k != 0
                    || !matches!(
                        linear_solve_policy,
                        BvpSciLinearSolvePolicy::AutoBanded
                            | BvpSciLinearSolvePolicy::ExperimentalBorderedBanded
                    )
                {
                    *calc_statistics
                        .entry("bvp sci direct banded route rejected".to_string())
                        .or_insert(0) += 1;
                    singular = true;
                    break;
                }

                let (df_dy, df_dp) = banded_jac_fn(x, &y, &p);
                if df_dp.is_some() {
                    *calc_statistics
                        .entry("bvp sci direct banded parameter derivative rejected".to_string())
                        .or_insert(0) += 1;
                    singular = true;
                    break;
                }
                let mut x_middle = faer_col::zeros(m - 1);
                for j in 0..(m - 1) {
                    x_middle[j] = x[j] + 0.5 * h[j];
                }
                let (df_dy_middle, df_dp_middle) = banded_jac_fn(&x_middle, &y_middle, &p);
                if df_dp_middle.is_some() {
                    *calc_statistics
                        .entry("bvp sci direct banded parameter derivative rejected".to_string())
                        .or_insert(0) += 1;
                    singular = true;
                    break;
                }
                let (dbc_dya, dbc_dyb, dbc_dp) = if let Some(bc_jac_fn) = bc_jac {
                    bc_jac_fn(&ya, &yb, &p)
                } else {
                    estimate_bc_jac(bc, &ya, &yb, &p, Some(&bc_res))
                };
                if dbc_dp.is_some() {
                    *calc_statistics
                        .entry("bvp sci direct banded boundary parameter rejected".to_string())
                        .or_insert(0) += 1;
                    singular = true;
                    break;
                }

                let assembly_started = Instant::now();
                let blocks = match assemble_bordered_banded_blocks_from_pointwise(
                    h,
                    &df_dy,
                    &df_dy_middle,
                    &dbc_dya,
                    &dbc_dyb,
                ) {
                    Ok(blocks) => blocks,
                    Err(_) => {
                        *calc_statistics
                            .entry("bvp sci direct banded assembly failures".to_string())
                            .or_insert(0) += 1;
                        singular = true;
                        break;
                    }
                };
                add_elapsed_micros_stat(
                    calc_statistics,
                    "bvp sci direct banded assembly us",
                    assembly_started,
                );
                *calc_statistics
                    .entry("bvp sci direct banded assembly calls".to_string())
                    .or_insert(0) += 1;
                *calc_statistics
                    .entry("bvp sci global sparse jacobian bypasses".to_string())
                    .or_insert(0) += 1;
                *calc_statistics
                    .entry("bvp sci linear route bordered banded".to_string())
                    .or_insert(0) += 1;

                custom_timer.jac_tac();
                custom_timer.linear_system_tic();
                let factor_started = Instant::now();
                let factorization = match factor_bordered_banded_structured(&blocks) {
                    Ok(factorization) => factorization,
                    Err(_) => {
                        *calc_statistics
                            .entry("bvp sci bordered factorization failures".to_string())
                            .or_insert(0) += 1;
                        singular = true;
                        break;
                    }
                };
                add_elapsed_micros_stat(
                    calc_statistics,
                    "bvp sci bordered factorization us",
                    factor_started,
                );
                *calc_statistics
                    .entry("bvp sci bordered factorization calls".to_string())
                    .or_insert(0) += 1;
                let mut workspace = factorization.new_workspace();
                let solve_started = Instant::now();
                let step_col = match factorization.solve_with_workspace(&res, &mut workspace) {
                    Ok(step_col) => step_col,
                    Err(_) => {
                        *calc_statistics
                            .entry("bvp sci bordered structured solve failures".to_string())
                            .or_insert(0) += 1;
                        singular = true;
                        break;
                    }
                };
                add_elapsed_micros_stat(
                    calc_statistics,
                    "bvp sci bordered structured solve us",
                    solve_started,
                );
                *calc_statistics
                    .entry("number of solving linear systems".to_string())
                    .or_insert(0) += 1;
                *calc_statistics
                    .entry("bvp sci linear backend bordered structured solves".to_string())
                    .or_insert(0) += 1;
                *calc_statistics
                    .entry("bvp sci bordered initial solve calls".to_string())
                    .or_insert(0) += 1;
                cost = step_col.squared_norm_l2();
                lu_solver = None;
                banded_solver = None;
                bordered_factorization = Some(factorization);
                bordered_workspace = Some(workspace);
            } else {
                // Compute sparse Jacobians
                let (df_dy, df_dp) = if let Some(jac_fn) = fun_jac {
                    jac_fn(x, &y, &p)
                } else {
                    estimate_fun_jac(fun, x, &y, &p, Some(&f))
                };

                let mut x_middle = faer_col::zeros(m - 1);
                for j in 0..(m - 1) {
                    x_middle[j] = x[j] + 0.5 * h[j];
                }

                let (df_dy_middle, df_dp_middle) = if let Some(jac_fn) = fun_jac {
                    jac_fn(&x_middle, &y_middle, &p)
                } else {
                    estimate_fun_jac(fun, &x_middle, &y_middle, &p, Some(&f_middle))
                };

                let (dbc_dya, dbc_dyb, dbc_dp) = if let Some(bc_jac_fn) = bc_jac {
                    bc_jac_fn(&ya, &yb, &p)
                } else {
                    estimate_bc_jac(bc, &ya, &yb, &p, Some(&bc_res))
                };

                // Construct global Jacobian (use cached path if available)
                let jac_matrix: faer_mat = if let Some(cache) = sparsity_cache.as_mut() {
                    construct_global_jac_cached(
                        n,
                        m,
                        k,
                        h,
                        &df_dy,
                        &df_dy_middle,
                        df_dp.as_ref().map(|v| v.as_slice()),
                        df_dp_middle.as_ref().map(|v| v.as_slice()),
                        &dbc_dya,
                        &dbc_dyb,
                        dbc_dp.as_ref(),
                        cache,
                    );
                    // Clone the matrix from the cache after values have been written
                    let cached_mat = cache.matrix_clone();

                    cached_mat
                } else {
                    construct_global_jac(
                        n,
                        m,
                        k,
                        h,
                        &df_dy,
                        &df_dy_middle,
                        df_dp.as_ref().map(|v| v.as_slice()),
                        df_dp_middle.as_ref().map(|v| v.as_slice()),
                        &dbc_dya,
                        &dbc_dyb,
                        dbc_dp.as_ref(),
                    )
                };
                custom_timer.jac_tac();
                custom_timer.linear_system_tic();
                let jacobian_diagnostics = inspect_global_jacobian(&jac_matrix);
                record_global_jacobian_diagnostics(calc_statistics, &jacobian_diagnostics);
                if jacobian_diagnostics.is_definitely_singular() {
                    *calc_statistics
                        .entry("global jacobian explicit singular precheck failures".to_string())
                        .or_insert(0) += 1;
                    singular = true;
                    break;
                }
                lu_solver = None;
                banded_solver = None;
                bordered_factorization = None;
                bordered_workspace = None;

                let route_profile = profile_bordered_banded_global_jacobian(
                    &jac_matrix,
                    n,
                    m,
                    k,
                    BvpSciBandedRoutePolicy::default(),
                );
                record_linear_backend_route(calc_statistics, &route_profile);
                let should_try_bordered = match linear_solve_policy {
                    BvpSciLinearSolvePolicy::ExperimentalBorderedBanded => {
                        if route_profile.recommended_route != BvpSciBandedRoute::BorderedBanded {
                            *calc_statistics
                                .entry("bvp sci bordered structured route rejected".to_string())
                                .or_insert(0) += 1;
                            singular = true;
                            break;
                        }
                        true
                    }
                    BvpSciLinearSolvePolicy::AutoBanded => {
                        let eligible = k == 0
                            && route_profile.recommended_route == BvpSciBandedRoute::BorderedBanded;
                        if k > 0 {
                            *calc_statistics
                                .entry("bvp sci auto bordered parameter fallback".to_string())
                                .or_insert(0) += 1;
                        }
                        eligible
                    }
                    _ => false,
                };
                let should_try_banded = match linear_solve_policy {
                    BvpSciLinearSolvePolicy::Sparse => false,
                    BvpSciLinearSolvePolicy::AutoBanded => {
                        route_profile.recommended_route == BvpSciBandedRoute::FullScalarBanded
                    }
                    BvpSciLinearSolvePolicy::RequireFullBanded => {
                        if route_profile.recommended_route != BvpSciBandedRoute::FullScalarBanded {
                            *calc_statistics
                                .entry("bvp sci full banded route rejected".to_string())
                                .or_insert(0) += 1;
                            singular = true;
                            break;
                        }
                        true
                    }
                    BvpSciLinearSolvePolicy::ExperimentalBorderedBanded => false,
                };
                let auto_banded_sparse_fallback = linear_solve_policy
                    == BvpSciLinearSolvePolicy::AutoBanded
                    && !should_try_banded
                    && !should_try_bordered;

                if should_try_bordered {
                    let extract_started = Instant::now();
                    let extracted_blocks = extract_bordered_banded_blocks(&jac_matrix, n, m, k);
                    add_elapsed_micros_stat(
                        calc_statistics,
                        "bvp sci bordered extraction us",
                        extract_started,
                    );
                    *calc_statistics
                        .entry("bvp sci bordered extraction calls".to_string())
                        .or_insert(0) += 1;
                    match extracted_blocks.and_then(|blocks| {
                        let factor_started = Instant::now();
                        let factorization = factor_bordered_banded_structured(&blocks);
                        add_elapsed_micros_stat(
                            calc_statistics,
                            "bvp sci bordered factorization us",
                            factor_started,
                        );
                        *calc_statistics
                            .entry("bvp sci bordered factorization calls".to_string())
                            .or_insert(0) += 1;
                        let factorization = factorization?;
                        let mut workspace = factorization.new_workspace();
                        let solve_started = Instant::now();
                        let step_col = factorization.solve_with_workspace(&res, &mut workspace);
                        add_elapsed_micros_stat(
                            calc_statistics,
                            "bvp sci bordered structured solve us",
                            solve_started,
                        );
                        *calc_statistics
                            .entry("bvp sci bordered structured solve calls".to_string())
                            .or_insert(0) += 1;
                        let step_col = step_col?;
                        Ok((factorization, workspace, step_col))
                    }) {
                        Ok((factorization, workspace, step_col)) => {
                            bordered_factorization = Some(factorization);
                            bordered_workspace = Some(workspace);
                            cost = step_col.squared_norm_l2();
                            *calc_statistics
                                .entry("number of solving linear systems".to_string())
                                .or_insert(0) += 1;
                            *calc_statistics
                                .entry(
                                    "bvp sci linear backend bordered structured solves".to_string(),
                                )
                                .or_insert(0) += 1;
                            *calc_statistics
                                .entry("bvp sci bordered initial solve calls".to_string())
                                .or_insert(0) += 1;
                        }
                        Err(_) if linear_solve_policy == BvpSciLinearSolvePolicy::AutoBanded => {
                            *calc_statistics
                                .entry("bvp sci bordered structured solve failures".to_string())
                                .or_insert(0) += 1;
                            *calc_statistics
                                .entry("bvp sci bordered factor fallback to sparse".to_string())
                                .or_insert(0) += 1;
                            match jac_matrix.sp_lu() {
                                Ok(lu) => {
                                    let step = lu.solve(res.as_mat());
                                    let step_col =
                                        faer_col::from_fn(step.nrows(), |i| *step.get(i, 0));
                                    lu_solver = Some(lu);
                                    cost = step_col.squared_norm_l2();
                                    *calc_statistics
                                        .entry("number of solving linear systems".to_string())
                                        .or_insert(0) += 1;
                                    *calc_statistics
                                        .entry(
                                            "bvp sci linear backend sparse fallback solves"
                                                .to_string(),
                                        )
                                        .or_insert(0) += 1;
                                }
                                Err(_) => {
                                    *calc_statistics
                                        .entry(
                                            "global jacobian lu factorization failures".to_string(),
                                        )
                                        .or_insert(0) += 1;
                                    singular = true;
                                    break;
                                }
                            }
                        }
                        Err(_) => {
                            *calc_statistics
                                .entry("bvp sci bordered structured solve failures".to_string())
                                .or_insert(0) += 1;
                            singular = true;
                            break;
                        }
                    }
                } else if should_try_banded {
                    match factor_banded_and_solve(&jac_matrix, &res) {
                        Ok((lu, step_col)) => {
                            banded_solver = Some(lu);
                            cost = step_col.squared_norm_l2();
                            *calc_statistics
                                .entry("number of solving linear systems".to_string())
                                .or_insert(0) += 1;
                            *calc_statistics
                                .entry("bvp sci linear backend full banded solves".to_string())
                                .or_insert(0) += 1;
                        }
                        Err(_) if linear_solve_policy == BvpSciLinearSolvePolicy::AutoBanded => {
                            *calc_statistics
                                .entry("bvp sci full banded factor fallback to sparse".to_string())
                                .or_insert(0) += 1;
                            match jac_matrix.sp_lu() {
                                Ok(lu) => {
                                    let step = lu.solve(res.as_mat());
                                    let step_col =
                                        faer_col::from_fn(step.nrows(), |i| *step.get(i, 0));

                                    lu_solver = Some(lu);
                                    cost = step_col.squared_norm_l2();
                                    *calc_statistics
                                        .entry("number of solving linear systems".to_string())
                                        .or_insert(0) += 1;
                                    *calc_statistics
                                        .entry(
                                            "bvp sci linear backend sparse fallback solves"
                                                .to_string(),
                                        )
                                        .or_insert(0) += 1;
                                }
                                Err(_) => {
                                    *calc_statistics
                                        .entry(
                                            "global jacobian lu factorization failures".to_string(),
                                        )
                                        .or_insert(0) += 1;
                                    singular = true;
                                    break;
                                }
                            }
                        }
                        Err(_) => {
                            *calc_statistics
                                .entry("bvp sci full banded factorization failures".to_string())
                                .or_insert(0) += 1;
                            singular = true;
                            break;
                        }
                    }
                } else {
                    // Attempt sparse LU decomposition
                    match jac_matrix.sp_lu() {
                        Ok(lu) => {
                            let step = lu.solve(res.as_mat());
                            let step_col = faer_col::from_fn(step.nrows(), |i| *step.get(i, 0));

                            lu_solver = Some(lu);
                            cost = step_col.squared_norm_l2();
                            *calc_statistics
                                .entry("number of solving linear systems".to_string())
                                .or_insert(0) += 1;
                            let counter_key = if auto_banded_sparse_fallback {
                                "bvp sci linear backend sparse fallback solves"
                            } else {
                                "bvp sci linear backend sparse solves"
                            };
                            *calc_statistics.entry(counter_key.to_string()).or_insert(0) += 1;
                        }
                        Err(_) => {
                            *calc_statistics
                                .entry("global jacobian lu factorization failures".to_string())
                                .or_insert(0) += 1;
                            singular = true;
                            break;
                        }
                    }
                }
            }
            njev += 1;
            custom_timer.linear_system_tac();
        }
        custom_timer.linear_system_tic();
        let step_col = if let Some(ref lu) = lu_solver {
            let step: Mat<f64> = lu.solve(res.as_mat());
            Some(faer_col::from_fn(step.nrows(), |i| *step.get(i, 0)))
        } else if let Some(ref lu) = banded_solver {
            match solve_with_banded_factor(lu, &res) {
                Ok(step) => Some(step),
                Err(_) => {
                    *calc_statistics
                        .entry("bvp sci full banded solve failures".to_string())
                        .or_insert(0) += 1;
                    singular = true;
                    None
                }
            }
        } else if let (Some(factorization), Some(workspace)) =
            (bordered_factorization.as_ref(), bordered_workspace.as_mut())
        {
            let solve_started = Instant::now();
            let step_result = factorization.solve_with_workspace(&res, workspace);
            add_elapsed_micros_stat(
                calc_statistics,
                "bvp sci bordered structured solve us",
                solve_started,
            );
            *calc_statistics
                .entry("bvp sci bordered structured solve calls".to_string())
                .or_insert(0) += 1;
            *calc_statistics
                .entry("bvp sci bordered reuse solve calls".to_string())
                .or_insert(0) += 1;
            match step_result {
                Ok(step) => Some(step),
                Err(_) => {
                    *calc_statistics
                        .entry("bvp sci bordered structured solve failures".to_string())
                        .or_insert(0) += 1;
                    singular = true;
                    None
                }
            }
        } else {
            None
        };
        custom_timer.linear_system_tac();

        if let Some(step_col) = step_col {
            // Extract step components
            let mut y_step = faer_dense_mat::zeros(n, m);
            for j in 0..m {
                for i in 0..n {
                    *y_step.get_mut(i, j) = step_col[j * n + i];
                }
            }

            let p_step = if k > 0 {
                faer_col::from_fn(k, |i| step_col[(m - 1) * n + (n + k) - k + i])
            } else {
                faer_col::zeros(0)
            };

            // Backtracking line search
            let mut alpha = 1.0;
            let mut best_y = y.clone();
            let mut best_p = p.clone();
            let mut best_cost = cost;

            info!("\n \n Starting backtracking line search \n \n");

            for trial in 0..=n_trial {
                if trial > 0 {
                    info!("\n \n Line search trial {}, alpha = {} \n \n", trial, alpha);
                }
                let mut y_new = y.clone();
                let mut p_new = p.clone();

                // y_new = y - alpha * y_step
                for i in 0..n {
                    for j in 0..m {
                        *y_new.get_mut(i, j) -= alpha * y_step.get(i, j);
                    }
                }

                // p_new = p - alpha * p_step
                for i in 0..k {
                    p_new[i] -= alpha * p_step[i];
                }
                custom_timer.fun_tic();
                // Compute new residuals
                let (col_res_new, _y_middle_new, _f_new, _f_middle_new) =
                    collocation_fun(fun, &y_new, &p_new, x, h);
                custom_timer.fun_tac();
                let ya_new = faer_col::from_fn(n, |i| *y_new.get(i, 0));
                let yb_new = faer_col::from_fn(n, |i| *y_new.get(i, m - 1));
                let bc_res_new = bc(&ya_new, &yb_new, &p_new);

                let mut res_new = faer_col::zeros((m - 1) * n + (n + k));
                for j in 0..(m - 1) {
                    for i in 0..n {
                        res_new[j * n + i] = *col_res_new.get(i, j);
                    }
                }
                for i in 0..(n + k) {
                    res_new[(m - 1) * n + i] = bc_res_new[i];
                }

                let step_new_col = if let Some(ref lu) = lu_solver {
                    let step_new = lu.solve(res_new.as_mat());
                    faer_col::from_fn(step_new.nrows(), |i| *step_new.get(i, 0))
                } else if let Some(ref lu) = banded_solver {
                    match solve_with_banded_factor(lu, &res_new) {
                        Ok(step) => step,
                        Err(_) => {
                            *calc_statistics
                                .entry("bvp sci full banded solve failures".to_string())
                                .or_insert(0) += 1;
                            singular = true;
                            break;
                        }
                    }
                } else if let (Some(factorization), Some(workspace)) =
                    (bordered_factorization.as_ref(), bordered_workspace.as_mut())
                {
                    let solve_started = Instant::now();
                    let step_result = factorization.solve_with_workspace(&res_new, workspace);
                    add_elapsed_micros_stat(
                        calc_statistics,
                        "bvp sci bordered structured solve us",
                        solve_started,
                    );
                    *calc_statistics
                        .entry("bvp sci bordered structured solve calls".to_string())
                        .or_insert(0) += 1;
                    *calc_statistics
                        .entry("bvp sci bordered line search solve calls".to_string())
                        .or_insert(0) += 1;
                    match step_result {
                        Ok(step) => step,
                        Err(_) => {
                            *calc_statistics
                                .entry("bvp sci bordered structured solve failures".to_string())
                                .or_insert(0) += 1;
                            singular = true;
                            break;
                        }
                    }
                } else {
                    singular = true;
                    break;
                };
                let cost_new = step_new_col.squared_norm_l2();

                if cost_new < (1.0 - 2.0 * alpha * sigma) * cost {
                    best_y = y_new;
                    best_p = p_new;
                    best_cost = cost_new;
                    break;
                }

                if trial < n_trial {
                    alpha *= tau;
                }
            }

            if singular {
                break;
            }

            y = best_y;
            p = best_p;
            cost = best_cost;

            // Check convergence
            let (col_res_final, _, _f_final, f_middle_final) = collocation_fun(fun, &y, &p, x, h);
            let ya_final = faer_col::from_fn(n, |i| *y.get(i, 0));
            let yb_final = faer_col::from_fn(n, |i| *y.get(i, m - 1));
            let bc_res_final = bc(&ya_final, &yb_final, &p);

            let mut converged = true;

            // Check collocation residuals
            for j in 0..(m - 1) {
                for i in 0..n {
                    if col_res_final.get(i, j).abs()
                        >= tol_r[j] * (1.0 + f_middle_final.get(i, j).abs())
                    {
                        converged = false;
                        break;
                    }
                }
                if !converged {
                    break;
                }
            }

            // Check boundary condition residuals
            for i in 0..(n + k) {
                if bc_res_final[i].abs() >= bc_tol {
                    converged = false;
                    break;
                }
            }

            if converged {
                info!("\n \n Newton iteration converged! \n \n");
                break;
            } else {
                info!("\n \n Newton iteration not converged, continuing... \n \n");
            }

            // Decide whether to recompute Jacobian
            recompute_jac = if alpha == 1.0 {
                false // Continue with same Jacobian
            } else {
                true // Recompute Jacobian
            };

            if njev >= max_njev {
                break;
            }
        } else {
            singular = true;
            break;
        }
    }

    (y, p, singular)
}

/// Create a cubic spline given values and derivatives at mesh nodes.
///
/// This function creates a cubic spline representation matching SciPy's PPoly format.
/// The formulas for coefficients are taken from scipy.interpolate.CubicSpline.
///
/// # Arguments
/// * `y` - Function values at mesh nodes (n, m)
/// * `yp` - Function derivatives at mesh nodes (n, m)
/// * `x` - Mesh nodes (m,)
/// * `h` - Mesh intervals (m-1,)
///
/// # Returns
/// PPoly spline with axis=1 for compatibility with SciPy
pub fn create_spline(y: &faer_dense_mat, yp: &faer_dense_mat, x: &faer_col, h: &faer_col) -> PPoly {
    let (n, m) = (y.nrows(), y.ncols());

    // Initialize coefficient array: c[degree][interval][polynomial]
    // For cubic spline: degree = 4 (coefficients for x^3, x^2, x^1, x^0)
    let mut c = vec![vec![vec![0.0; n]; m - 1]; 4];

    // Compute cubic spline coefficients for each interval and each polynomial
    for j in 0..(m - 1) {
        for i in 0..n {
            // Compute slope over interval
            let slope = (y.get(i, j + 1) - y.get(i, j)) / h[j];

            // Compute t parameter: (yp_left + yp_right - 2*slope) / h
            let t = (yp.get(i, j) + yp.get(i, j + 1) - 2.0 * slope) / h[j];

            // Cubic spline coefficients (highest degree first)
            c[0][j][i] = t / h[j]; // x^3 coefficient
            c[1][j][i] = (slope - yp.get(i, j)) / h[j] - t; // x^2 coefficient
            c[2][j][i] = *yp.get(i, j); // x^1 coefficient
            c[3][j][i] = *y.get(i, j); // x^0 coefficient
        }
    }

    PPoly {
        c,
        x: (0..x.nrows()).map(|i| x[i]).collect(),
        extrapolate: Extrapolate::Bool(true),
        axis: 1, // axis=1 for compatibility with SciPy
    }
}
/// Estimate RMS values of collocation residuals using Lobatto quadrature.
///
/// The residuals are defined as the difference between the derivatives of
/// our solution and RHS of the ODE system. We use relative residuals, i.e.,
/// normalized by 1 + abs(f). RMS values are computed as sqrt from the
/// normalized integrals of the squared relative residuals over each interval.
/// Integrals are estimated using 5-point Lobatto quadrature.
///
/// # Arguments
/// * `fun` - ODE function f(x, y, p)
/// * `sol` - Solution spline (PPoly)
/// * `x` - Mesh nodes (m,)
/// * `h` - Mesh intervals (m-1,)
/// * `p` - Parameters (k,)
/// * `r_middle` - Residuals at middle points (n, m-1)
/// * `f_middle` - Function values at middle points (n, m-1)
///
/// # Returns
/// RMS residuals for each interval (m-1,)
pub fn estimate_rms_residuals(
    fun: &ODEFunction,
    sol: &PPoly,
    x: &faer_col,
    h: &faer_col,
    p: &faer_col,
    r_middle: &faer_dense_mat,
    f_middle: &faer_dense_mat,
) -> faer_col {
    let (n, m1) = (r_middle.nrows(), r_middle.ncols()); // n equations, m-1 intervals
    info!("\n \n Estimating RMS residuals using Lobatto quadrature \n \n");

    // Compute middle points of intervals
    let mut x_middle = faer_col::zeros(m1);
    for j in 0..m1 {
        x_middle[j] = x[j] + 0.5 * h[j];
    }

    // Compute Lobatto quadrature points: x_middle ± s
    // where s = 0.5 * h * sqrt(3/7)
    let mut x1 = faer_col::zeros(m1);
    let mut x2 = faer_col::zeros(m1);
    for j in 0..m1 {
        let s = 0.5 * h[j] * ((3.0 / 7.0) as f64).sqrt();
        x1[j] = x_middle[j] + s;
        x2[j] = x_middle[j] - s;
    }

    // Evaluate solution and its derivative at quadrature points
    let x1_vec: Vec<f64> = (0..x1.nrows()).map(|i| x1[i]).collect();
    let x2_vec: Vec<f64> = (0..x2.nrows()).map(|i| x2[i]).collect();

    let y1 = sol.call(&x1_vec, &[m1], Some(0), None);
    let y2 = sol.call(&x2_vec, &[m1], Some(0), None);
    let y1_prime = sol.call(&x1_vec, &[m1], Some(1), None);
    let y2_prime = sol.call(&x2_vec, &[m1], Some(1), None);

    // Convert to faer matrices
    let y1_faer = faer_dense_mat::from_fn(y1.ncols(), y1.nrows(), |i, j| y1[(j, i)]);
    let y2_faer = faer_dense_mat::from_fn(y2.ncols(), y2.nrows(), |i, j| y2[(j, i)]);

    // Evaluate ODE RHS at quadrature points
    let f1 = fun(&x1, &y1_faer, p);
    let f2 = fun(&x2, &y2_faer, p);

    // Compute residuals: r = y' - f
    let mut r1 = faer_dense_mat::zeros(n, m1);
    let mut r2 = faer_dense_mat::zeros(n, m1);
    for j in 0..m1 {
        for i in 0..n {
            *r1.get_mut(i, j) = y1_prime[(j, i)] - f1.get(i, j);
            *r2.get_mut(i, j) = y2_prime[(j, i)] - f2.get(i, j);
        }
    }

    // Normalize residuals by (1 + |f|) and compute squared norms
    let mut r_middle_norm = faer_col::zeros(m1);
    let mut r1_norm = faer_col::zeros(m1);
    let mut r2_norm = faer_col::zeros(m1);

    for j in 0..m1 {
        let mut sum_r_middle = 0.0;
        let mut sum_r1 = 0.0;
        let mut sum_r2 = 0.0;

        for i in 0..n {
            // Normalize residuals
            let r_mid_norm = r_middle.get(i, j) / (1.0 + f_middle.get(i, j).abs());
            let r1_norm_val = r1.get(i, j) / (1.0 + f1.get(i, j).abs());
            let r2_norm_val = r2.get(i, j) / (1.0 + f2.get(i, j).abs());

            // Sum of squares (for complex numbers, this would be |r|^2)
            sum_r_middle += r_mid_norm * r_mid_norm;
            sum_r1 += r1_norm_val * r1_norm_val;
            sum_r2 += r2_norm_val * r2_norm_val;
        }

        r_middle_norm[j] = sum_r_middle;
        r1_norm[j] = sum_r1;
        r2_norm[j] = sum_r2;
    }

    // Apply 5-point Lobatto quadrature formula and take square root
    let mut rms_res = faer_col::zeros(m1);
    for j in 0..m1 {
        let integral =
            0.5 * (32.0 / 45.0 * r_middle_norm[j] + 49.0 / 90.0 * (r1_norm[j] + r2_norm[j]));
        rms_res[j] = integral.sqrt();
    }

    rms_res
}
/// Modify mesh by inserting nodes
///
/// # Arguments
/// * `x` - Current mesh nodes
/// * `insert_1` - Intervals to insert 1 node in middle
/// * `insert_2` - Intervals to insert 2 nodes (divide into 3 parts)
///
/// # Returns
/// New mesh with inserted nodes
pub fn modify_mesh(x: &faer_col, insert_1: &[usize], insert_2: &[usize]) -> faer_col {
    let mut new_points: Vec<f64> = (0..x.nrows()).map(|i| x[i]).collect();

    // Insert 1 node in middle of intervals
    for &i in insert_1 {
        let mid = 0.5 * (x[i] + x[i + 1]);
        new_points.push(mid);
    }

    // Insert 2 nodes to divide interval into 3 parts
    for &i in insert_2 {
        let p1 = (2.0 * x[i] + x[i + 1]) / 3.0;
        let p2 = (x[i] + 2.0 * x[i + 1]) / 3.0;
        new_points.push(p1);
        new_points.push(p2);
    }

    new_points.sort_by(|a, b| a.partial_cmp(b).unwrap());
    faer_col::from_fn(new_points.len(), |i| new_points[i])
}

/// BVP solver with mesh refinement.
///
/// This compatibility wrapper uses the default Newton/line-search strategy.
pub fn solve_bvp(
    fun: &ODEFunction,
    bc: &BCFunction,
    x: faer_col,
    y: faer_dense_mat,
    p: Option<faer_col>,
    _s: Option<faer_dense_mat>,
    fun_jac: Option<&ODEJacobian>,
    bc_jac: Option<&BCJacobian>,
    tol: f64,
    max_nodes: usize,
    verbose: u8,
    bc_tol: Option<f64>,
    customtimer: Option<CustomTimer>,
) -> Result<BVPResult, String> {
    solve_bvp_with_strategy_params(
        fun,
        bc,
        x,
        y,
        p,
        _s,
        fun_jac,
        bc_jac,
        tol,
        max_nodes,
        verbose,
        bc_tol,
        customtimer,
        None,
    )
}

/// BVP solver with mesh refinement and explicit Newton/line-search strategy.
pub fn solve_bvp_with_strategy_params(
    fun: &ODEFunction,
    bc: &BCFunction,
    x: faer_col,
    y: faer_dense_mat,
    p: Option<faer_col>,
    _s: Option<faer_dense_mat>,
    fun_jac: Option<&ODEJacobian>,
    bc_jac: Option<&BCJacobian>,
    tol: f64,
    max_nodes: usize,
    verbose: u8,
    bc_tol: Option<f64>,
    customtimer: Option<CustomTimer>,
    strategy_params: Option<&HashMap<String, Vec<f64>>>,
) -> Result<BVPResult, String> {
    solve_bvp_with_strategy_and_linear_policy(
        fun,
        bc,
        x,
        y,
        p,
        _s,
        fun_jac,
        bc_jac,
        tol,
        max_nodes,
        verbose,
        bc_tol,
        customtimer,
        strategy_params,
        BvpSciLinearSolvePolicy::Sparse,
    )
}

/// BVP solver with explicit Newton strategy and linear backend policy.
pub fn solve_bvp_with_strategy_and_linear_policy(
    fun: &ODEFunction,
    bc: &BCFunction,
    x: faer_col,
    y: faer_dense_mat,
    p: Option<faer_col>,
    _s: Option<faer_dense_mat>,
    fun_jac: Option<&ODEJacobian>,
    bc_jac: Option<&BCJacobian>,
    tol: f64,
    max_nodes: usize,
    verbose: u8,
    bc_tol: Option<f64>,
    customtimer: Option<CustomTimer>,
    strategy_params: Option<&HashMap<String, Vec<f64>>>,
    linear_solve_policy: BvpSciLinearSolvePolicy,
) -> Result<BVPResult, String> {
    solve_bvp_with_strategy_linear_policy_and_banded_jacobian(
        fun,
        bc,
        x,
        y,
        p,
        _s,
        fun_jac,
        None,
        bc_jac,
        tol,
        max_nodes,
        verbose,
        bc_tol,
        customtimer,
        strategy_params,
        linear_solve_policy,
    )
}

/// BVP solver entry point with an optional native pointwise banded Jacobian.
///
/// Supplying `banded_fun_jac` selects the direct bordered-banded assembly route
/// for parameter-free endpoint problems and bypasses global sparse Jacobian
/// materialization.
pub fn solve_bvp_with_strategy_linear_policy_and_banded_jacobian(
    fun: &ODEFunction,
    bc: &BCFunction,
    mut x: faer_col,
    mut y: faer_dense_mat,
    p: Option<faer_col>,
    _s: Option<faer_dense_mat>,
    fun_jac: Option<&ODEJacobian>,
    banded_fun_jac: Option<&ODEBandedJacobian>,
    bc_jac: Option<&BCJacobian>,
    tol: f64,
    max_nodes: usize,
    verbose: u8,
    bc_tol: Option<f64>,
    customtimer: Option<CustomTimer>,
    strategy_params: Option<&HashMap<String, Vec<f64>>>,
    linear_solve_policy: BvpSciLinearSolvePolicy,
) -> Result<BVPResult, String> {
    if fun_jac.is_none() {
        info!("\n \n No Jacobian provided, using numerical estimation \n \n");
    } else {
        info!("\n \n Analytical Jacobian provided! \n \n");
    }

    //  info!("\n \n Starting BVP solver with {} equations and {} mesh points \n \n", n, m);
    let mut custom_timer = customtimer.unwrap_or_else(|| CustomTimer::new());
    let vec_of_tuples = vec![
        ("number of iterations".to_string(), 0),
        ("number of solving linear systems".to_string(), 0),
        ("number of jacobians recalculations".to_string(), 0),
        ("number of grid refinements".to_string(), 0),
    ];
    let mut calc_statistics: HashMap<String, usize> = vec_of_tuples.into_iter().collect();

    let n = y.nrows();
    let mut m = x.nrows();
    let singular_plumbing = match _s.as_ref() {
        Some(singular_term) => Some(SingularTermPlumbing::new(singular_term, x[0])?),
        None => None,
    };

    if y.ncols() != m {
        return Err("y must have same number of columns as x has elements".to_string());
    }

    let mut p = p.unwrap_or_else(|| faer_col::zeros(0));
    let bc_tol = bc_tol.unwrap_or(tol);
    let max_iteration = 10;

    if let Some(plumbing) = singular_plumbing.as_ref() {
        plumbing.project_initial_guess(&mut y);
        *calc_statistics
            .entry("bvp sci singular term enabled".to_string())
            .or_insert(0) += 1;
    }

    let singular_enabled = singular_plumbing.is_some();
    let fun_wrapped: Option<Box<ODEFunction>> = singular_plumbing.as_ref().map(|plumbing| {
        let plumbing = plumbing.clone();
        let fun_ptr = fun as *const ODEFunction;
        Box::new(
            move |x: &faer_col, y: &faer_dense_mat, p: &faer_col| -> faer_dense_mat {
                let base = unsafe { (&*fun_ptr)(x, y, p) };
                plumbing.wrap_rhs_output(x, &base, y)
            },
        ) as Box<ODEFunction>
    });
    let fun_ref: &ODEFunction = fun_wrapped
        .as_ref()
        .map(|wrapped| wrapped.as_ref())
        .unwrap_or(fun);
    let fun_jac_ref = if singular_enabled { None } else { fun_jac };
    let banded_fun_jac_ref: Option<&ODEBandedJacobian> = if singular_plumbing.is_some() {
        None
    } else {
        banded_fun_jac
    };

    custom_timer.fun_tic();
    // Initial validation
    let f_test = fun_ref(&x, &y, &p);
    if (f_test.nrows(), f_test.ncols()) != (y.nrows(), y.ncols()) {
        return Err(format!(
            "Function return shape ({}, {}) doesn't match y shape ({}, {})",
            f_test.nrows(),
            f_test.ncols(),
            y.nrows(),
            y.ncols()
        ));
    }
    custom_timer.fun_tac();
    let ya_test = faer_col::from_fn(n, |i| *y.get(i, 0));
    let yb_test = faer_col::from_fn(n, |i| *y.get(i, m - 1));
    let bc_test = bc(&ya_test, &yb_test, &p);
    let expected_bc_size = n + p.nrows();
    if bc_test.nrows() != expected_bc_size {
        return Err(format!(
            "BC return size {} doesn't match expected size {}",
            bc_test.nrows(),
            expected_bc_size
        ));
    }

    let mut status = 0;
    let mut iteration = 0;

    if verbose == 2 {
        info!(
            "{:^15}{:^15}{:^15}{:^15}{:^15}",
            "Iteration", "Max residual", "Max BC residual", "Total nodes", "Nodes added"
        );
    }

    // Create sparsity cache for Jacobian reuse across Newton iterations
    let k = p.nrows();
    let mut sparsity_cache: Option<JacobianSparsityCache> = None;

    loop {
        m = x.nrows();

        // Compute mesh intervals
        let mut h = faer_col::zeros(m - 1);
        for i in 0..(m - 1) {
            h[i] = x[i + 1] - x[i];
        }

        // Recreate cache if mesh size changed (m changes after grid refinement)
        let needs_new_cache = match &sparsity_cache {
            Some(cache) => cache.m != m || cache.n != n || cache.k != k,
            None => true,
        };
        if needs_new_cache {
            sparsity_cache = Some(JacobianSparsityCache::new(n, m, k));
        }

        // Solve Newton system
        let (y_new, p_new, singular) = solve_newton(
            n,
            m,
            &h,
            fun_ref,
            bc,
            fun_jac_ref,
            banded_fun_jac_ref,
            bc_jac,
            y.clone(),
            p.clone(),
            &x,
            tol,
            bc_tol,
            &mut custom_timer,
            &mut calc_statistics,
            sparsity_cache.as_mut(),
            strategy_params,
            linear_solve_policy,
        );

        y = y_new;
        p = p_new;
        *calc_statistics
            .entry("number of iterations".to_string())
            .or_insert(0) += 1;
        iteration += 1;

        // Compute collocation residuals and boundary condition residuals
        let (col_res, _y_middle, f, f_middle) = collocation_fun(fun_ref, &y, &p, &x, &h);
        let ya_curr = faer_col::from_fn(n, |i| *y.get(i, 0));
        let yb_curr = faer_col::from_fn(n, |i| *y.get(i, m - 1));
        let bc_res = bc(&ya_curr, &yb_curr, &p);
        let max_bc_res = (0..bc_res.nrows())
            .map(|i| bc_res[i].abs())
            .fold(0.0, f64::max);

        if singular {
            status = 2;
            break;
        }

        // This relation is not trivial, but can be verified
        let mut r_middle = faer_dense_mat::zeros(n, m - 1);
        for j in 0..(m - 1) {
            for i in 0..n {
                *r_middle.get_mut(i, j) = 1.5 * col_res.get(i, j) / h[j];
            }
        }

        let sol = create_spline(&y, &f, &x, &h);
        custom_timer.fun_tic();
        let rms_res = estimate_rms_residuals(fun, &sol, &x, &h, &p, &r_middle, &f_middle);
        custom_timer.fun_tac();
        let max_rms_res = (0..rms_res.nrows()).map(|i| rms_res[i]).fold(0.0, f64::max);

        // Determine which intervals need refinement
        let mut insert_1 = Vec::new();
        let mut insert_2 = Vec::new();

        for j in 0..(m - 1) {
            if rms_res[j] > tol && rms_res[j] < 100.0 * tol {
                insert_1.push(j);
            } else if rms_res[j] >= 100.0 * tol {
                insert_2.push(j);
            }
        }

        let nodes_added = insert_1.len() + 2 * insert_2.len();

        if m + nodes_added > max_nodes {
            status = 1;
            if verbose == 2 {
                info!(
                    "{:^15}{:^15.2e}{:^15.2e}{:^15}{:^15}",
                    iteration,
                    max_rms_res,
                    max_bc_res,
                    m,
                    format!("({})", nodes_added)
                );
            }
            break;
        }

        if verbose == 2 {
            info!(
                "{:^15}{:^15.2e}{:^15.2e}{:^15}{:^15}",
                iteration, max_rms_res, max_bc_res, m, nodes_added
            );
        }

        if nodes_added > 0 {
            custom_timer.grid_refinement_tic();
            x = modify_mesh(&x, &insert_1, &insert_2);
            // Evaluate solution at new mesh points
            let x_eval: Vec<f64> = (0..x.nrows()).map(|i| x[i]).collect();
            let y_new_vals = sol.call(&x_eval, &[x.nrows()], Some(0), None);
            y = faer_dense_mat::from_fn(y_new_vals.ncols(), y_new_vals.nrows(), |i, j| {
                y_new_vals[(j, i)]
            });
            if let Some(plumbing) = singular_plumbing.as_ref() {
                plumbing.project_initial_guess(&mut y);
            }
            custom_timer.grid_refinement_tac();
            *calc_statistics
                .entry("number of grid refinements".to_string())
                .or_insert(0) += 1;
        } else if max_bc_res <= bc_tol {
            status = 0;
            break;
        } else if iteration >= max_iteration {
            status = 3;
            break;
        }
    }

    if verbose > 0 {
        match status {
            0 => info!(
                "Solved in {} iterations, number of nodes {}.",
                iteration,
                x.nrows()
            ),
            1 => log::warn!("Number of nodes exceeded after iteration {}.", iteration),
            2 => log::warn!("Singular Jacobian encountered on iteration {}.", iteration),
            3 => log::warn!(
                "Unable to satisfy boundary conditions tolerance on iteration {}.",
                iteration
            ),
            _ => {}
        }
    }
    custom_timer.fun_tic();
    let final_f = fun_ref(&x, &y, &p);
    let final_h = faer_col::from_fn(x.nrows() - 1, |i| x[i + 1] - x[i]);
    let (col_res_final, _y_middle_final, _, f_middle_final) =
        collocation_fun(fun_ref, &y, &p, &x, &final_h);
    custom_timer.fun_tac();
    let mut r_middle_final = faer_dense_mat::zeros(n, x.nrows() - 1);
    for j in 0..(x.nrows() - 1) {
        for i in 0..n {
            *r_middle_final.get_mut(i, j) = 1.5 * col_res_final.get(i, j) / final_h[j];
        }
    }
    let final_sol = create_spline(&y, &final_f, &x, &final_h);
    custom_timer.fun_tic();
    let final_rms_res = estimate_rms_residuals(
        fun,
        &final_sol,
        &x,
        &final_h,
        &p,
        &r_middle_final,
        &f_middle_final,
    );
    custom_timer.fun_tac();
    let message = match status {
        0 => "The algorithm converged to the desired accuracy.",
        1 => "The maximum number of mesh nodes is exceeded.",
        2 => "A singular Jacobian encountered when solving the collocation system.",
        3 => "The solver was unable to satisfy boundary conditions tolerance on iteration 10.",
        _ => "Unknown status",
    };

    Ok(BVPResult {
        sol: Some(final_sol),
        p: if p.nrows() > 0 { Some(p) } else { None },
        x,
        y,
        yp: final_f,
        rms_residuals: final_rms_res,
        niter: iteration,
        status,
        message: message.to_string(),
        success: status == 0,
        calc_statistics: calc_statistics,
        custom_timer: custom_timer,
        strategy_params: Some(
            strategy_params
                .cloned()
                .unwrap_or_else(default_strategy_params),
        ),
    })
}

/// Compatibility alias for the production faer sparse BVP solver.
///
/// Historically this function duplicated the full solve loop. Keep the public
/// name for callers, but route through the single implementation so strategy,
/// timings, mesh refinement, and Jacobian-cache behavior cannot diverge.
pub fn solve_bvp_sparse(
    fun: &ODEFunction,
    bc: &BCFunction,
    x: faer_col,
    y: faer_dense_mat,
    p: Option<faer_col>,
    _s: Option<faer_dense_mat>,
    fun_jac: Option<&ODEJacobian>,
    bc_jac: Option<&BCJacobian>,
    tol: f64,
    max_nodes: usize,
    verbose: u8,
    bc_tol: Option<f64>,
    customtimer: Option<CustomTimer>,
) -> Result<BVPResult, String> {
    solve_bvp_with_strategy_params(
        fun,
        bc,
        x,
        y,
        p,
        _s,
        fun_jac,
        bc_jac,
        tol,
        max_nodes,
        verbose,
        bc_tol,
        customtimer,
        None,
    )
}