I have a problem in my code. Who can help me to fix it?.. I can't find my mistake. Please help me.

close; clear; clc;
lambda = 532e-9;  % laser Nd:Yag wavelength second harmonic
k = 2 .*pi ./lambda;
w0 = 0.01;    % Laser radius
z = 5000; %.* (0.01:0.005:1);  % 5 km distance
[x,y] = meshgrid((-1:0.01:1).*0.5);
[phi,r] = cart2pol(x,y);    % Changing the polar coordinate variable to Cartesian
z_R = k .*w0 .^2 ./2;    % Constant
w = w0 .*sqrt(1 + (z./z_R) .^2);    % Constant
R = z_R .*(z./z_R + z_R./z);    % Constant
l = 1;    % l = ..., -3, -2, -1, 0, 1, 2, 3, ...
p = 2;    % p = 0, 1, 2, 3, ...
N = abs(l) + 2 .*p;
Clp = sqrt(2 .*factorial(p) ./(pi .*factorial(p + abs(l))));   % clp is a parameter in laguerre gaussian beam equation that you can see in " https://en.wikipedia.org/wiki/Gaussian_beam "
% This below laguerre0-3 is lagurre polynomials that you can see in " https://en.wikipedia.org/wiki/Laguerre_polynomials "
laguerre0 = '@(A,l) 1';
laguerre1 = '@(A,l) -A + l + 1';
laguerre2 = '@(A,l) A.^2 ./2 - (l + 2) .*A + (l + 1) .*(l + 2) ./2';
laguerre3 = '@(A,l) -1 .*A.^3 ./6 + (l + 3) .*A.^2 ./2 - (l + 2) .*(l +3) .*A ./2 + (l + 1) .*(l + 2) .*(l+3) ./6';
fun = str2func("laguerre" + num2str(p));
% u is laguerre-gaussian equation you can see in " https://en.wikipedia.org/wiki/Gaussian_beam "
u = Clp ./w .*(sqrt(2 .*(x.^2 + y.^2)) ./w) .^abs(l) .*exp(-1 .*(x.^2 + y.^2) ./w .^2) ...
    .*fun(2 .*(x.^2 + y.^2) ./w.^2,l) .*exp(-1i .*k .*(x.^2 + y.^2) ./(2 .*R)) ...
    .*exp(1i .*l .*phi) .*exp(1i .*(N + 1) .*atan(z ./z_R));
I = abs(u).^2;
s=pcolor(x,y,I);
s.EdgeColor = "none";

Answer 1

[Stephen23]

Rather than forcing pseudo-indices into variable names and then attempting to use STR2FUNC.... simply use a cell array:

lambda = 532e-9;  % laser Nd:Yag wavelength second harmonic
k = 2 .*pi ./lambda;
w0 = 0.01;    % Laser radius
z = 5000; %.* (0.01:0.005:1);  % 5 km distance
[x,y] = meshgrid((-1:0.01:1).*0.5);
[phi,r] = cart2pol(x,y);    % Changing the polar coordinate variable to Cartesian
z_R = k .*w0 .^2 ./2;    % Constant
w = w0 .*sqrt(1 + (z./z_R) .^2);    % Constant
R = z_R .*(z./z_R + z_R./z);    % Constant
l = 1;    % l = ..., -3, -2, -1, 0, 1, 2, 3, ...
p = 2;    % p = 0, 1, 2, 3, ...
N = abs(l) + 2 .*p;
Clp = sqrt(2 .*factorial(p) ./(pi .*factorial(p + abs(l))));   % clp is a parameter in laguerre gaussian beam equation that you can see in " https://en.wikipedia.org/wiki/Gaussian_beam "
% This below laguerre0-3 is lagurre polynomials that you can see in " https://en.wikipedia.org/wiki/Laguerre_polynomials "
laguerre = {...
    @(A,l) 1,...
    @(A,l) -A + l + 1,...
    @(A,l) A.^2 ./2 - (l + 2) .*A + (l + 1) .*(l + 2) ./2,...
    @(A,l) -1 .*A.^3 ./6 + (l + 3) .*A.^2 ./2 - (l + 2) .*(l +3) .*A ./2 + (l + 1) .*(l + 2) .*(l+3) ./6};
fun = laguerre{p+1};
% u is laguerre-gaussian equation you can see in " https://en.wikipedia.org/wiki/Gaussian_beam "
u = Clp ./w .*(sqrt(2 .*(x.^2 + y.^2)) ./w) .^abs(l) .*exp(-1 .*(x.^2 + y.^2) ./w .^2) ...
    .*fun(2 .*(x.^2 + y.^2) ./w.^2,l) .*exp(-1i .*k .*(x.^2 + y.^2) ./(2 .*R)) ...
    .*exp(1i .*l .*phi) .*exp(1i .*(N + 1) .*atan(z ./z_R));
I = abs(u).^2;
s = pcolor(x,y,I);
s.EdgeColor = "none";

[Babr]

Thank you so much.

[Stephen23]

@Babr: please do not forget to click the accept button if this answer helped you!

[Babr]

just, in line 20 " fun = laguerre{p+1} "

This way it becomes more complete.

[Stephen23]

@Babr: good point! I corrected my answer.