Double integral of a function (Beginners Question)

Since there are a lot of floats in the integrand, I'm guessing
an approximate result is OK. Instead of

> integrate(integrand(kappa,Cn2,(1-xi)/L,L,Lambda,Theta_inv,xi,k,r),kappa,0,inf);
> integrate(%,xi,0,1);
> float(%);

try this:

I : integrand(kappa,Cn2,(1-xi)/L,L,Lambda,Theta_inv,xi,k,r);
I1 : quad_qags (I, xi, 0, 1)[1];
I2 : quad_qagi (I1, kappa, 0, inf);

 => [3.93347674966401e-20, 9.416575845421898e-22, 3225, 5]

The integrand is very nearly zero almost everywhere,
which might make it more difficult for the quadrature functions
to estimate the error. Maybe you can rescale the integrand.

Note that I2 doesn't take into account any error estimate from I1.
Iterated quadrature is probably not such a great method,
although Monte Carlo and quasi-Monte Carlo are probably less
efficient (i.e. need more function evaluations for the same
accuracy) than iterated quadrature for a 2-d problem.
Maybe there are better methods especially for low-dimensional
integration.

plot3d (I, [kappa, 0.1, 1], [xi, 0, 1], [grid, 50, 50]);
seems to be informative.

By the way, what is the origin of this problem?

best

Robert Dodier