I’d like to calculate the following double integral
Integrate[exp(-x)Integrate[exp(-y),{y,0,x}],{x,0,∞}]
I try the following in julia but i think i probably get it wrong.
using QuadGK
inner_function(y) = exp(-y)
outer_function(x) = quadgk(y → inner_function(y), 0, x)[1]
result, _ = quadgk(outer_function, 0, Inf)
println(“Integration result:”, result)
Would anyone help me with it? Thanks!
Hi! 
First off, why do you think you got it wrong? It’s very useful to add what kind of errors or unexpected results you already found and it makes it easier for other people to help.
Second, I think the outer_function is missing multiplication by exp(-x) like so
outer_function(x) = exp(-x) * quadgk(inner_function, 0, x)[1]
Now this is the whole integrand you wrote down (note that you can shorten y -> f(y) to just f in general).
Other than that, it looks correct to me. After the change, I get
julia> result, _ = quadgk(outer_function, 0, Inf)
(0.49999999999999994, 4.8229187783801426e-11)
thank you! this is so nice! I didn’t realize i was almost there.
Also note that for performance (if it matters), you generally want to set a tolerance for both integrals and use a stricter tolerance for the inner integral by a factor of ~10 to make sure that the outer integral doesn’t try to shape the inter integral.