In MATLAB, *quadgk(@(s) (exp(-s.^2)-1)./s.^2,0,1)* gives me 0.8615, consistent with the correct answer, 1 - 1/e -sqrt(pi) erf(1).

In Julia, *quadgk(s-> exp(-s.^2-1.0)./s.^2,0.0,1.0)* gives

me

*ERROR: DomainError with 4.7733380679681323e-153:*

*integrand produced Inf in the interval (0.0, 9.546676135936265e-153)*

which suggests that it’s getting very hung up about the behaviour near s=0, even though Gaussian quadrature isn’t meant to evaluate the integrand at s=0.

Of course, this isn’t the actual integral that I’m interested in, but it illustrates the point.

Am I missing something subtle here? I really want to try using Julia, but this is not a great start.