How to plot a numerical integration?

Hello guys, I made a numerical integration of this function using adaptive Gauss-Kronrod quadrature 'quadgk()':

as you can see the integrate is k, but for the purpose that I want, I need to plot that result vs z, which is the depth. So, for that I made this code:

    %%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.
    z = linspace(0,z1,100000);

    %%numerical integration

    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 = quadgk(fun,0,Inf);

    %%Plot of the total fluence rate, collimated and diffuse fluence rate

    PHIc = exp(-U1t.*z);
    PHI = PHIc + PHId;
    figure;
    semilogy(z,PHI)
    xlabel('Reduced depth z/l*');
    ylabel('Normalized Fluence Rate \phi');

Basically, I made a vector from 0 to z1 (which is the value that I need to reach) and once the numerical integration is made I tried to plot it vs z. I do not know if I am making something wrong, I am new in MATLAB and I am really thankful for your answers.

Answer 1

[Aquatris]

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');

[Johan Sebastian Diaz Tovar]

Thank you so much for your soon answer!! Can I ask u something else? Due that I'm new in MATLAB, I would like to know how to make a nested for, it is because I would like to plot, something like this:

And following your idea it will be making a nested for, for the variable r. Thank u one more time!

[Aquatris]

I'll assume the r on the plot is r0 in the code. Then what you need to do is;

for i = 1:length(zSet)
   for j = 1:length(rSet)
      z  = zSet(i);
      r0 = rSet(j);
      % rest of the code 
      %
      %
      PHIc(i,j) = exp(-U1t.*z);
      PHI(i,j) = PHIc(i,j) + PHId(i,j);
   end
end

Then you need to use meshgrid function;

[Z,R] = meshgrid(zSet,rSet); % creates X-Y grid
contourf(Z,R,PHI) % creates plot

Sometimes I confuse the order of Z,R so you might wanna check if x-axis of the plot is really zSet values!

I also recommend you to follow the examples from the contour() (contourf() is a filled version of contour()) function documentation, which can be found here, to better understand the function and its use.

[Johan Sebastian Diaz Tovar]

Thank you so much!!!! I am really thankful with your answers!

[Aquatris]

Happy to help. Hope these are what you were looking for.