The problem is, the quadgk function expect a scalar value from your function. However you assign z to be a vector, hence your function is not scalar. I modified your code and reduce the step size for z to reduce the computational time to be able to test it. Here is my modifications for your code, feel free to ask anything if you do not understand any part.
%%Total fluence parameters and constants
nu = 0; %zeroth-order Bessel subscript for the function of the first kind.
g = 0.9; %anisotropy coefficient
f = g.^2; g1 = g./(g+1); %first two moments of the delta P1 phase function
n = 1.4; A = -0.13755*n.^3 + 4.3390*n.^2 -4.90366*n + 1.6896; %refraction index and cubic polynomial estimation
%for FRC
Ups = 0.060; Ua = 0.0576; Us = Ups./(1-g); %Reduced, absorption and scattering coefficients
Utr = Ups + Ua; h = 2./(3*Utr);
Ueff = sqrt(3*Ua*Utr); U1s = Us*(1 - f);
U1t = Ua + U1s; z1 = Utr./Ups; %z is the direction of the collimated light within the medium.
r0 = 3./Utr; %Gaussian beam radius.
zSet = linspace(0,z1,100);
%%numerical integration
for i = 1:length(zSet)
z = zSet(i);
fun = @(x) (((3.*U1s.*(U1t + g.*Ua).*r0.^2.*exp(-(r0.*x).^2)./8)./4.*(x.^2 + Ueff.^2 - U1t.^2)).*exp(-Utr.*z) + ...
((-3.*g1.*U1s.*r0.^2.*exp(-(r0.*x).^2)./8 - 4.*(1./(A.*h) + U1t.^2).* ...
((3.*U1s.*(U1t + g.*Ua).*r0.^2.*exp(-(r0.*x).^2)./8)./4.*(x.^2 + Ueff.^2 - U1t.^2)))./4.*(1./(A.*h) + sqrt(x.^2 + ...
Ueff.^2))).* ...
exp(-sqrt(x.^2+Ueff.^2).*z)).*besselj(nu,x).*x;
PHId(i) = quadgk(fun,0,Inf);
%%Plot of the total fluence rate, collimated and diffuse fluence rate
PHIc(i) = exp(-U1t.*z);
PHI(i) = PHIc(i) + PHId(i);
end
figure;
semilogy(zSet,PHI)
xlabel('Reduced depth z/l*');
ylabel('Normalized Fluence Rate \phi');
