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.
