.png)
"Split and conquer" is always a good idea for a complex problem like this. When trying to do everything at once, the chances of getting it right are rather small. Split into smaller parts and test to make sure they make sense before moving on. Here is a possible start:
Instead of using a double integral, try for instance to first calculate the B function:
function B = attiliB(alpha,x)
B = zeros(size(x));
for i = 1:length(x)
fun = @(y) sqrt(1-y.^2).*besselj(0,alpha*x(i).*y).*y;
B(i) = integral(fun,0,1);
end
Test it by:
>> b = 25.2;
>> h = 50;
>> x = linspace(0,1);
>> z = attiliB(a,x);
>> plot(x,z)
Now you can use attiliB inside a single integral for the I_ksi function. Doing this instead of a double integral may perhaps be slower, but the chances of getting it right are certainly better.
It turns out that B(x) in this case is very well approximated by a second order polynomial, which you can use as a slightly less accurate, but certainly much faster, substitute for attiliB:
>> pB = polyfit(x,z,2);
>> plot(x,z,x,polyval(pB,x),'.')
Good luck with the rest!
