The problem is that norminv(x,0,1) goes to -Inf as x goes to zero and Inf as x goes to 1 (and ditto for y and z), so it's hard to integrate accurately. If you can manage with a larger tolerance,
q3 = integral3(fun3,0,1,0,1,0,1,'AbsTol',1e-3);
does not return any warnings (the default tolerance is 1e-10). If you need higher precision, you'll just need to be patient, as the warnings are not fatal. However, I don't know if the precision is actually met under such circumstances.
