ApproxFun does this:
julia> f = Fun(x->exp(-x^2),0..5);
julia> g = 2/sqrt(π)*integrate(f); h = g - g(0);
julia> x = 0:0.5:5; [f.(x) h.(x)]
11×2 Matrix{Float64}:
1.0 0.0
0.778801 0.5205
0.367879 0.842701
0.105399 0.966105
0.0183156 0.995322
0.00193045 0.999593
0.00012341 0.999978
4.78512e-6 0.999999
1.12535e-7 1.0
1.60523e-9 1.0
1.3888e-11 1.0
For a univariate function, I doubt you will be able to do much better.