Is it possible to curve fit with a custom model that contains an integration of a Bessel function?

Is it possible to use the following equation for "B" as a custom model in the curve fitting toolbox? The only fitting parameter is "ilam".

T=9;
Tc=9.25;
z0=1;
xi0=0.04;
xiv=xi0/sqrt(1-(T/Tc));
lam=0.04;
ilam=1/lam; % ilam is going to be my fitting parameter.
phi0=20.7;
x0=-5:0.1:5;
%
for i=1:length(x0)
    x=x0(i);
    scale=1;
    B(i)=(phi0/(2*pi))*quad((@(q) (besselk(1,(sqrt(q.^2+ilam^2).*xiv),scale)/((sqrt(q.^2+ilam^2)+q).*lam.*besselk(1,(xiv/lam),scale))).*exp(-q.*z0).*besselj(0,q.*x).*q),0,20);
end
%
figure
hold all
plot(x0,B,'b')
xlabel('X (\mum)')
ylabel('B (G)')
box on