Function compute_diagonal_jacobian 

pub fn compute_diagonal_jacobian<F, Func>(
    f: &Func,
    t: F,
    y: &Array1<F>,
    f_current: &Array1<F>,
    block_size: usize,
) -> Array2<F>
where
    F: IntegrateFloat,
    Func: Fn(F, ArrayView1<'_, F>) -> Array1<F>,

Computes Jacobian for a system with diagonal or block-diagonal structure