MATLAB going from Cart to Pol back to Cart coords for cylindrical plot

1. What I WANT to do:
(i) Use input n to generate an n*n cartesian grid

[x y] = meshgrid(linspace(-1,1,n));  

(ii) Generate polar coordinates

[theta r] = cart2pol(x,y);  

(iii) Evaluate a function in cylindrical coordinates

z = f(theta,r);  

(iv) Plot the result using (say) pcolor (or surf, or anything)

pcolor(x,y,abs(z).^2) %The function is complex, a Laguerre-Gauss to be exact.  

2. What I CAN do... The only way I can get the plots to work is by starting with my polar parameters and working back to cartesian from there:
(i) Define parameters

r=linspace(0,1,n); theta=linspace(0,2*pi,n);  

(ii) Create both grids and evaluate f

[theta r]=meshgrid(theta,r);  
[x y]=pol2cart(theta,r);  
z=f(theta,r);  

(iii) Plot

pcolor(x,y,abs(z).^2)  

The PROBLEM is that now my grid is circular, and I would like to evaluate the function everywhere ON A RECTANGULAR grid (because my analysis depends on having square pixel arrays). The reiterate, using method 2 above, I get a circular plot circumscribed in a square; imagine a black circle with white along the edges... but I WANT to evaluate the function in this "white" region. HOWEVER, using method 1 does NOT work -- the function is all messed up when I plot (Just google Laguerre-Gauss modes to see what the plots should look like).

I want to be able to start with a rect grid and assign every point a polar coordinate, instead of start with polar coordinates and assign them all cartesian points.

I've been messing with this on an off for a long time, and I can't figure out how to get around this seemingly simple issue.

Edit 1
It seems that the problem lies in how the coordinate matrices are generated. Below I've posted screen-shots of a simple 3by3 example illustrating how approach 1 and approach 2 generate different numbers.

How to make these numbers compatible?

I have no reputation points so I cannot upload the images directly... links below show the 3by3 example... see comments for links to actual images of the Laguerre-Gauss plots I'm trying to make...

apply cart2pol
apply pol2cart

Edit 2

Currently, the result of approach (1.) gives wrong results, as shown here:

desired approach, wrong result

The second approach gives the right images, unfortunately it's only a circle and not the entire square. It is shown here:

implemented approach, limited result

3D plots of both approaches are shown here - only the colorful part of the top figure is correct.

Edit 3

Here is a screenshot of the function f which is being used above. Note, that it asks for more input parameters than just r,theta. Typical values are:

w0 = 0.5;
p = 0;
l = 5;

The function C gives a normalization and L are Laguerre polynomials. Both of these functions have been thoroughly tested and yield the expected results.

Edit 4
Here is enough code to run my example z=U(0,5,r,phi,w0)+U(0,-5,r,phi,w0); explicitly. The plot itself is given by pcolor(x,y,abs(z).^2).

Note that the Lpl() function is inserted as a comment. This will have to be saved as its own m-file for the U function to run properly.

%% Laguerre-Gauss Modes U = U(p,l,r,phi,w0)
%  Source: OAM theory paper section 2.A eqn 1.
%  Assuming POLAR coordinates and evaluating AT beam waist.
% -- That is, z=0 for w(z)=w0(sqrt(1+z/zR)) 
% ---- ie, w(0) = w0
% Assuming z=0 also renders the Gouy phase arctan(z/zR) irrelevant.
% Note: Rayleigh Range zR is not explicitly defined because z=0 --> it is irrelevant too.
% Since zR is the only wavelength dependent term, wavelength also doesn't
% matter.

function out = U(p,l,r,phi,w0)
%Function handles for clarity
e = @(x) exp(x);
C = @(p,l) sqrt((2*factorial(p))/(pi*factorial(p+abs(l))));
L = @(p,l,z) Lpl(p,l,z);

%% Lpl() FUNCTION
% function out = Lpl(p,l,z)
% 
% l=abs(l);
% LL=0;
% for mm=1:p+1
%     m=mm-1;
%     L=LL;
%     LL= L+((-1)^m)*(factorial(p+l)/(factorial(p-m)*factorial(l+m)*factorial(m)))*(z.^m); 
% end
% out = LL;

%% 

out = (C(p,l)/w0)*...
    (((sqrt(2).*r)/w0)^abs(l))*...
    (e((-r.^2)/w0^2))*...
    (L(p,l,((2.*r.^2)/w0^2)))*...
    (e((-1)*1i*l.*phi)); ``

Edit
The answer was rewritten based on the code provided in Edit 4 of the question.

I think the trouble stems from the function U. You don't apply element wise operations to all parts of the equation. If you change it to:

out = (C(p,l)./w0).* ...              % here it's a .* instead of *
    (((sqrt(2).*r)./w0).^abs(l)).* ... % here it's a .* instead of *
    (e((-r.^2)./w0.^2)).* ...          % here it's a .* instead of *
    (L(p,l,((2.*r.^2)./w0.^2))).* ...  % here it's a .* instead of *
    (e((-1)*1i*l.*phi));

You get the following two results, shown below.

This figure used an input of cartesian coordinates:

And this figure used the polar coordinates:

The "coarser" resolution in the second figure is due to the less suitable resolution of the grid. But in essence you resolve the same features.