how can i improve time ( less time ) in matlab integration

how can i spend less time during the calculation ? I_value ( triple integral ) calculation

function I_value = doit
T=0.119;
s=0.319;
R=0.531;
tic
f2 = @(r,b,g) 1./(r.^2.*sqrt(1 - (b./r).^2 - (g^-2)*((2/15)*(s/(R))^9 *(1./(r - 1).^9 - 1./(r + 1).^9 - 9./(8* r).*(1./(r - 1).^8 - 1./(r + 1).^8)) - (s/(R))^3 *(1./(r -1).^3 - 1./(r + 1).^3 - 3./(2* r).* (1./(r - 1).^2 - 1./(r + 1).^2)))));
% The folloing function only works sith scalar b and g values.
X_scalar_b_scalar_g = @(b,g)real(pi - 2*b*quadgk(@(r)f2(r,b,g),rmin(g,b,R(k)),Inf,'AbsTol',1e-3,'RelTol',1e-3,'MaxIntervalCount',3000));
% Make X work with array inputs for b and a scalar g value.
X_scalar_g = @(b,g)arrayfun(@(b)X_scalar_b_scalar_g(b,g),b); 
f3 = @(b,g) 2*(1 - cos(X_scalar_g(b,g))).*b; 
qQd_scalar_g = @(g)quadgk(@(b)f3(b,g),0.001,15,'AbsTol',1e-3,'RelTol',1e-3,'MaxIntervalCount',7000);
% Make qQd_scalar_g work with array g inputs.
qQd = @(g)arrayfun(qQd_scalar_g,g);
f4 =@(g) g.^5.*qQd(g)./(exp(g.^2/T(i)));
I_value = (1/T^3)*quadgk(f4,0.07,2,'AbsTol',1e-6,'RelTol',1e-6,'MaxIntervalCount',7000)
toc
end
function r = rmin(b,g,R)
f1 =  @(r)  1 - (b./r).^2 - (g^-2)*((2/15)*(0.319/R)^9 *(1./(r - 1).^9 - 1./(r + 1).^9 - 9./(8*r).*(1./(r - 1).^8 - 1./(r + 1).^8)) -(0.319/R)^3 *(1./(r-1).^3 - 1./(r+1).^3 - 3./(2*r).*(1./(r-1).^2 - 1./(r+1).^2)));
r = fzero(f1,[1.00000001,50]);
end

It takes a lot of minutes. Can I change it to C, via some program ? Ideas ?