How to calculate double integral?

Hi! I have a problem with solution of double integral. The syms solves too long and it cannot be used.But in another way I have problem.

I have an integral fun2 on z and it has x which is a variable in the integral fun0. How can I set a variable x to first calculate the integral f2 over z, and then integral f3 over x? (Of course, when I set x=some number I obtain a curve or a set of curves if x=0:0.1:1 and make for j = 1:length(x), but I doubt about this result, because the behavior of curves is not correct).

clear all, close all

n=1;

t=1;

r=1;

s=0:0.2:10;

for i = 1:length(s)

k=s(i);

fun2=@(z)(z.*exp(2.*n.*t.*z.^2).*(besselj(0,(k.*z.*x)))./sqrt(1-z.^2));

f2(i,:)=integral(fun2,0,1);

fun0=@(x)(((((x.^2.*exp(-2.*t.*x.^2)./(x.^2+1/(r.^2)).^2)))).*f2(i));

f3(i,:)=integral(fun0,0,inf);

end

Cor=8/(r*(pi)^(3/2))*sqrt(2*n*t)*exp(-2*n*t)/(erf(sqrt(2*n*t))*((1+4*t/r^2)*exp(2*t/r^2)*erfc(sqrt(2*t/r^2))-2*sqrt(2*t)/(r*sqrt(pi)))).*f3;

plot(s,Cor,'b-');

Answer 1

[Torsten]

n = 1 ;
t = 1;
r = 1;
s = 0:0.2:10;
fun = @(x,z,k) x.^2.*exp(-2.*t.*x.^2)./(x.^2+1/r^2).^2 .* z.*exp(2*n*t.*z.^2).*besselj(0,k.*z.*x)./sqrt(1-z.^2);
f3 = arrayfun(@(k)integral2(@(x,z)fun(x,z,k),0,Inf,0,1),s)

f3 = 1×51

0.2959 0.2948 0.2915 0.2860 0.2784 0.2690 0.2579 0.2454 0.2316 0.2169 0.2015 0.1857 0.1697 0.1538 0.1381 0.1229 0.1083 0.0945 0.0815 0.0694 0.0583 0.0483 0.0392 0.0311 0.0240 0.0177 0.0123 0.0077 0.0038 0.0005

Cor = 8/(r*(pi)^(3/2))*sqrt(2*n*t)*exp(-2*n*t)/(erf(sqrt(2*n*t))*((1+4*t/r^2)*exp(2*t/r^2)*erfc(sqrt(2*t/r^2))-2*sqrt(2*t)/(r*sqrt(pi))))*f3

Cor = 1×51

1.0000 0.9962 0.9849 0.9663 0.9409 0.9091 0.8716 0.8292 0.7828 0.7331 0.6811 0.6276 0.5735 0.5196 0.4667 0.4153 0.3659 0.3192 0.2753 0.2346 0.1972 0.1631 0.1324 0.1051 0.0809 0.0598 0.0416 0.0259 0.0127 0.0017

plot(s,Cor,'b-')
grid on

[Hexe]

Dear Torsten, thank you very much :) It solves my problem and now I will know what to do with integrals like these.