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";
