Error computing a 3d integral

I am trying to compute a 3d integral by first computing a 2d integral using integral2 and then integral for the outer integral. I get an error :

Matrix diemnsions must agree

when I run the code below. The 2d integral is computed correctly by vzf.

Thanks !

---------------

rp=0.5;

h = 0.19;

zp = h;

huc =0.0;

vzj = vzf(0.7,0.5,0.5,0.0)

dv = integral(@(r) vzf(r,h,h,huc)*r,0,1.0)

function velzz =vzf(rp,zp,h,huc)

%function vzf is the 2d integral

velzz = integral2(@(rho,theta) velzi(rho,theta,rp,zp,h,huc),huc,1.0,0,2*pi);

end

function fz = velzi(rho,theta,rp,zp,h,huc)

% integrand of the 2d integral

fz = pload(rho,huc).*rho.*(zp+h)./(rp.^2 + rho.^2 + (zp+h).^2 + 2*rho.*rp.*cos(theta)).^(3/2) ;

end

function p = pload(r,huc)

p = 1.0*(r-huc).*(2 - r - huc);

end

[darova]

Do you have a source formulas? Image or something?