Integral problem numerical error

Hi, I am trying to solve a very complex equation numerically. I have following problem. Does anyone solve my problem? Thanks in advance.. My codes and error are in following.

function Fig1_pc_Gauss_scint_vs_norm_source_size clc

global kpp; global ksi; global s1x; global s1y; global s2x; global s2y;

PrT=7; %Prandtl numbers for temperature PrS=700; %Prandtl numbers for salinity beta=0.72; %Obukhov-Corrsin constant Q=2.75;

dist=100; %distance

KT=1; %Thermal diffusivity KS=1; %salt diffusivity XT=1e-8; %dissipation rate of mean squared temperature lambda=0.55e-6; %wavelength w=-3; %relative strength of temperature and salinity eta=1e-3; %kolmogorov microscale length eps=1e-7;

alphas=0.03;

px=0; py=0; %receiver coordinate

teta=KT/KS;

A=2.6e-4; B=1.75e-4;

k=2*pi/lambda;

j=1;

for rx=0:0.001:0.05

ry=rx;

r=(rx*rx+ry*ry)^0.5;

rt(j)=r*100;

coeff1=-1*pi*k*k*dist*beta*A*A*XT*(w-1)*(w-1)*((eps)^(-1/3))/(w*w*(w*w*teta+1-w*(1+teta))); cf1=1/(lambda*lambda*dist*dist);

dr1=@(s1x,s2x,s1y,s2y,kpp,ksi)((w*w*teta*(((exp(-3*beta*(eta^(4/3))/(2*PrT)*(kpp.^(4/3))-(Q*beta*eta*eta/PrT)*kpp.*kpp)).*(kpp.^(-8/3))))+... (((exp(-3*beta*(eta^(4/3))*(PrT+PrS)/(4*PrT*PrS)*(kpp.^(4/3))-(Q*beta*eta*eta*(PrT+PrS)/(2*PrT*PrS))*kpp.*kpp)).*(kpp.^(-8/3))))-... w*(1+teta)*(((exp(-3*beta*(eta^(4/3))/(2*PrS)*(kpp.^(4/3))-(Q*beta*eta*eta/PrS)*kpp.*kpp)).*(kpp.^(-8/3))))+... w*w*teta*Q*(eta^(2/3))*(((exp(-3*beta*(eta^(4/3))/(2*PrT)*(kpp.^(4/3))-(Q*beta*eta*eta/PrT)*kpp.*kpp)).*(kpp.^(-2))))+... Q*(eta^(2/3))*(((exp(-3*beta*(eta^(4/3))/(2*PrS)*(kpp.^(4/3))-(Q*beta*eta*eta/PrS)*kpp.*kpp)).*(kpp.^(-2))))-... w*(1+teta)*Q*(eta^(2/3))*(((exp(-3*beta*(eta^(4/3))*(PrT+PrS)/(4*PrT*PrS)*(kpp.^(4/3))-(Q*beta*eta*eta*(PrT+PrS)/(2*PrT*PrS))*kpp.*kpp)).*(kpp.^(-2))))).*... ((1-besselj(0,(kpp.*(((ksi*px+(1-ksi).*s1x-(1-ksi).*s2x).*(ksi*px+(1-ksi).*s1x-(1-ksi).*s2x)+(ksi*py+(1-ksi).*s1y-(1-ksi).*s2y).*(ksi*py+(1-ksi).*s1y-(1-ksi).*s2y)).^0.5))))));

innerIntegral1=@(s1x,s1y,s2x,s2y)(integral2(@(kpp,ksi)dr1(s1x,s1y,s2x,s2y,kpp,ksi),0,inf,0,1));

fc=@(s1x,s1y,s2x,s2y)((exp((-s1x.*s1x-s1y.*s1y)/(2*alphas*alphas))).*(exp((-s2x.*s2x-s2y.*s2y)/(2*alphas*alphas))).*... (exp((1i)*pi*((s1x-px).*(s1x-px)+(s1y-py).*(s1y-py)-(s2x-(px+rx)).*(s2x-(px+rx))-(s2y-(py+ry)).*(s2y-(py+ry)))/(lambda*dist))).*... (exp(coeff1*(innerIntegral1))));

innerIntegral2 = @(s1x)(@(y)integral3(@(s1y,s2x,s2y)fc(y,s1y,s2x,s2y),-inf,inf,-inf,inf,-inf,inf));

fcr(j) = (integral(innerIntegral2,-inf,inf))*cf1

j=j+1;

end

hold on plot(rt,fcr1,'--k','LineWidth',2); hold off

xlabel( '\it\bfr\rm (cm)','FontSize',24,'FontSize',24,'FontName','Times New Roman'); set(gca,'FontSize',24,'FontName','Times New Roman'); ylabel('|<\itu\rm(\bfp\rm)\itu*\rm(\bfp\rm+\bfr\rm)>|','FontSize',24,'FontName','Times New Roman');set(gca,'FontSize',24,'FontName','Times New Roman');

Error using integralCalc/finalInputChecks (line 511) Input function must return 'double' or 'single' values. Found 'function_handle'.

Error in integralCalc/iterateScalarValued (line 315) finalInputChecks(x,fx);

Error in integralCalc/vadapt (line 133) [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 104) [q,errbnd] = vadapt(@minusInfToInfInvTransform,interval);

Error in integral (line 89) Q = integralCalc(fun,a,b,opstruct);

Error in figure01_r_eksen_dongu_L_v1 (line 72) fcr(j) = (integral(innerIntegral2,-inf,inf))*cf1