It seems like it is lacking the functionality R's integrate has for handling infinite boundaries, but I suppose you could implement that yourself on the outside.

For what it's worth,

    use integrate::adaptive_quadrature::simpson::adaptive_simpson_method;
    use statrs::distribution::{Continuous, Normal};

    fn dnorm(x: f64) -> f64 {
        Normal::new(0.0, 1.0).unwrap().pdf(x)
    }
    
    fn main() {
        let result = adaptive_simpson_method(dnorm, -100.0, 100.0, 1e-2, 1e-8);
        println!("Result: {:?}", result);
    }

It does seem to be a usability hazard that the function being integrated is defined as a fn, rather than a Fn, as you can't pass closures that capture variables, requiring the weird dnorm definition