How to plot combined surf and 2d plot

I need to plot the intensity of a beam at z=1m as shown in figure.

I wrote code to find the intensity at z=1m. But failed to plot this curve. Code attaching below.

clc;clear all;close all;
x0 = -.03:0.00015:.03;
y0 = x0;
[x,y] = meshgrid(x0) ;
z=1;
lambda=532*10^(-9);
wo = 0.01;
k = (2.*pi)./(lambda);
Pho=0.0750  ;
Gamma = ((1i.*k)./(2.*z)) + (1./(Pho.^2)) + (1./(wo.^2));
Gamma_star = subs(Gamma, 1i, -1i);
alpha = (Gamma_star) - 1./((Gamma.*Pho.^4));
beta = 1 + 1./((Pho.^2).*Gamma);
beta_s = -1 + 1./((Pho.^2).*Gamma);
C = (1./(4.*(z.^2))).*((k.^2)./(alpha.*Gamma)) ;
fx = ((1i.*k)./(2.*z)) + ((k.*k)./(4.*z.*z.*Gamma)) +  ((k.*k.*beta.*beta)./(4.*z.*z.*alpha));
fy = ((-1i.*k)./(2.*z)) + ((k.*k)./(4.*z.*z.*Gamma)) +  ((k.*k.*beta.*beta)./(4.*z.*z.*alpha));
X= C.*exp( -(fx.*x.^2) - (fy.*y.^2) ) ;
I= abs(double((vpa(X))))./max(abs(double((vpa(X)))));
% to plot curve
Z = peaks(x,y) ;
C = I;
surfc(x,y,Z,C)
colorbar

[dpb]

The pertinent code reduces to...

x0 = -.03:0.00015:.03;
y0 = x0;
[x,y] = meshgrid(x0) ;
....
Z = peaks(x,y) ;
surfc(x,y,Z)

Over such a small range, the surface is essentially just a plane as is shown by the resulting plot -- you certainly can't expect the full range of the peaks function to be reproduced in just one tiny little portion of the area; it isn't a Mandelbrot set...

hold on
[X,Y,Z1]=peaks(25);                                 % get the default ranges used by peaks()
surfc(X,Y,Z1,'facecolor','none','linestyle',':')    % plot all transparently with only dashed lines 
zoom(2)                                             % blow up so can see original area

Now you can see your little part plotted in persepective with the whole thing (zoomed by 2X, so it's really only one-fourth, half in each direction centered around the origin), but you see you've just computed and drawn a very tiny piece out of the whole.

What all the other code has to do with anything or its purpose is totally unclear, but what you see is what you asked for as far as the plot goes...

[Athira T Das]

@dpb I tried to replot, but still have some error. I am trying to plot the beam profile at z=100m after propagating through a medium. I_numerical is the avergae intensity of the beam at z=100m, so I need to plot I_numerical with x,y.

clc;clear all;close all;
x=linspace(-8,8,500);
y=x;
[x,y] = meshgrid(x);
m =1;Omega =.1;

z=100 ; %propagation distance
sgm =1;wo = 0.01;
lambda=532 *10^-9;k=2*pi/lambda;
Pho=0.0162;Gamma = ((1i.*k)./(2.*z)) + (1./(Pho.^2)) + (1./(wo.^2));
Gamma_star = subs(Gamma, 1i, -1i);alpha = (Gamma_star) - 1./((Gamma.*Pho.^4));
beta = 1 + 1./((Pho.^2).*Gamma);beta_s = -1 + 1./((Pho.^2).*Gamma);
C = (1./(4.*(z.^2))).*((k.^2)./(alpha.*Gamma)) ;
fx = ((1i.*k)./(2.*z)) + ((k.*k)./(4.*z.*z.*Gamma)) +  ((k.*k.*beta.*beta)./(4.*z.*z.*alpha));
fy = ((-1i.*k)./(2.*z)) + ((k.*k)./(4.*z.*z.*Gamma)) +  ((k.*k.*beta.*beta)./(4.*z.*z.*alpha));
X= C.*exp( -(fx.*x.^2) - (fy.*y.^2) ) ;
I = (vpa(X));
I_numerical = abs(double(I))./max(abs(double(I))); % intensity at z=100m

surfc(x,y,I_numerical)
colorbar

●

[dpb]

x=linspace(-8,8,5);
y=x;
[x,y] = meshgrid(x);
m =1;
Omega=0.1;
z=100 ; %propagation distance
sgm =1;
wo = 0.01;
lambda=532 *10^-9;
k=2*pi/lambda;
Pho=0.0162;
Gamma = ((1i.*k)./(2.*z)) + (1./(Pho.^2)) + (1./(wo.^2));
Gamma_star = ((1i.*k)./(2.*z)) + (1./(Pho.^2)) + (1./(wo.^2));
%Gamma_star = subs(Gamma, 1i, -1i);
alpha = (Gamma_star) - 1./((Gamma.*Pho.^4));
beta  = 1 + 1./((Pho.^2).*Gamma);
beta_s=-1 + 1./((Pho.^2).*Gamma);
C = (1./(4.*(z.^2))).*((k.^2)./(alpha.*Gamma)) ;
fx=(( 1i.*k)./(2.*z)) + ((k.*k)./(4.*z.*z.*Gamma)) +  ((k.*k.*beta.*beta)./(4.*z.*z.*alpha));
fy=((-1i.*k)./(2.*z)) + ((k.*k)./(4.*z.*z.*Gamma)) +  ((k.*k.*beta.*beta)./(4.*z.*z.*alpha));
fx,fy

fx = 1.9442e+04 - 5.5760e+04i

fy = 1.9442e+04 - 1.7387e+05i

(fx.*x.^2),(fy.*y.^2) 

ans = 
   1.0e+06 *

   1.2443 - 3.5686i   0.3111 - 0.8922i   0.0000 + 0.0000i   0.3111 - 0.8922i   1.2443 - 3.5686i
   1.2443 - 3.5686i   0.3111 - 0.8922i   0.0000 + 0.0000i   0.3111 - 0.8922i   1.2443 - 3.5686i
   1.2443 - 3.5686i   0.3111 - 0.8922i   0.0000 + 0.0000i   0.3111 - 0.8922i   1.2443 - 3.5686i
   1.2443 - 3.5686i   0.3111 - 0.8922i   0.0000 + 0.0000i   0.3111 - 0.8922i   1.2443 - 3.5686i
   1.2443 - 3.5686i   0.3111 - 0.8922i   0.0000 + 0.0000i   0.3111 - 0.8922i   1.2443 - 3.5686i

ans = 
   1.0e+07 *

   0.1244 - 1.1127i   0.1244 - 1.1127i   0.1244 - 1.1127i   0.1244 - 1.1127i   0.1244 - 1.1127i
   0.0311 - 0.2782i   0.0311 - 0.2782i   0.0311 - 0.2782i   0.0311 - 0.2782i   0.0311 - 0.2782i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0311 - 0.2782i   0.0311 - 0.2782i   0.0311 - 0.2782i   0.0311 - 0.2782i   0.0311 - 0.2782i
   0.1244 - 1.1127i   0.1244 - 1.1127i   0.1244 - 1.1127i   0.1244 - 1.1127i   0.1244 - 1.1127i

X= C.*exp( -(fx.*x.^2) - (fy.*y.^2) );
V=abs(X)./max(abs(X));
I = (vpa(X));
I_numerical = abs(double(I))./max(abs(double(I))); % intensity at z=100m

%surfc(x,y,I_numerical)
%surfc(x,y,V)
V, I, I_numerical

V = 5×5

NaN NaN 0 NaN NaN NaN NaN 0 NaN NaN NaN NaN 1 NaN NaN NaN NaN 0 NaN NaN NaN NaN 0 NaN NaN

<mw-icon icon-id="meatballMenuUI" icon-width="16" icon-height="16" icon-config="{}" svg-path=""></mw-icon>

<mw-icon icon-id="kebabMenuUI" icon-width="16" icon-height="16" icon-config="{}" svg-path=""></mw-icon>

●

I = 

I_numerical = 5×5

NaN NaN 0 NaN NaN NaN NaN 0 NaN NaN NaN NaN 1 NaN NaN NaN NaN 0 NaN NaN NaN NaN 0 NaN NaN

<mw-icon icon-id="meatballMenuUI" icon-width="16" icon-height="16" icon-config="{}" svg-path=""></mw-icon>

<mw-icon icon-id="kebabMenuUI" icon-width="16" icon-height="16" icon-config="{}" svg-path=""></mw-icon>

%colorbar

We have no way to know what you're trying to model, but it appears the formulation is such that everything is NaN outside the one column of for which x=0 and the only finite value non-zero is for x,y=0.

Ensure you've actually written the equations correctly first, then examine the magnitudes of the exponential terms; as shown above, the exponential arguments being of magnitudes 10^6 and 10^7 is going to ensure they're going to underflow any reasonable numeric precision.

[Athira T Das]

@dpb The value is infinity if the range of x and y is from +8 to -8. Instead of that if the beam is plotted from -.03 to +0.03, the graph is perfect. The beam is defined only these ranges.

clc;clear all;close all;
x0 = -0.03:0.001:0.03;
y0 = x0;
[x,y] = meshgrid(x0) ;
z=1;
lambda=532*10^(-9);
wo = 0.01;
k = (2.*pi)./(lambda);
Pho=0.0750  ;
Gamma = ((1i.*k)./(2.*z)) + (1./(Pho.^2)) + (1./(wo.^2));
Gamma_star = subs(Gamma, 1i, -1i);
alpha = (Gamma_star) - 1./((Gamma.*Pho.^4));
beta = 1 + 1./((Pho.^2).*Gamma);
beta_s = -1 + 1./((Pho.^2).*Gamma);
C = (1./(4.*(z.^2))).*((k.^2)./(alpha.*Gamma)) ;
fx = ((1i.*k)./(2.*z)) + ((k.*k)./(4.*z.*z.*Gamma)) +  ((k.*k.*beta.*beta)./(4.*z.*z.*alpha));
fy = ((-1i.*k)./(2.*z)) + ((k.*k)./(4.*z.*z.*Gamma)) +  ((k.*k.*beta.*beta)./(4.*z.*z.*alpha));
X= C.*exp( -(fx.*x.^2) - (fy.*y.^2) ) ;
I= abs(double((vpa(X))))./max(abs(double((vpa(X)))));
plot(x0,I)

●

% to plot curve
Z = peaks(x,y) ;
C = I;
surfc(x,y,Z,C)
colorbar

[dpb]

clc;clear all;close all;
x0 = -0.03:0.001:0.03;
y0 = x0;
[x,y] = meshgrid(x0) ;
z=1;
lambda=532*10^(-9);
wo = 0.01;
k = (2.*pi)./(lambda);
Pho=0.0750  ;
Gamma = ((1i.*k)./(2.*z)) + (1./(Pho.^2)) + (1./(wo.^2));
Gamma_star = subs(Gamma, 1i, -1i);
alpha = (Gamma_star) - 1./((Gamma.*Pho.^4));
beta = 1 + 1./((Pho.^2).*Gamma);
beta_s = -1 + 1./((Pho.^2).*Gamma);
C = (1./(4.*(z.^2))).*((k.^2)./(alpha.*Gamma)) ;
fx = ((1i.*k)./(2.*z)) + ((k.*k)./(4.*z.*z.*Gamma)) +  ((k.*k.*beta.*beta)./(4.*z.*z.*alpha));
fy = ((-1i.*k)./(2.*z)) + ((k.*k)./(4.*z.*z.*Gamma)) +  ((k.*k.*beta.*beta)./(4.*z.*z.*alpha));
X= C.*exp( -(fx.*x.^2) - (fy.*y.^2) ) ;
I= abs(double((vpa(X))))./max(abs(double((vpa(X)))));
plot(x0,I)

% to plot curve
%Z = peaks(x,y) ;
%C = I;
surfc(x,y,I)
colorbar

●

Why do you keep bringing in the builtin function peaks? It has nothing whatsoever to do with anything you're trying to do other than being an example function used in the documentation to illustrate the 3D plotting functions and/or meshgrid.

The problem is that when you try to evaluate over such a large range outside the range of validity of the numerics, you aren't getting any points in the specific region of interest and so everything other than the origin is returned as NaN.

Note that if you set the limits of the plot so large, the ROI is so small in comparison to the range of the axes you've set that the result again looks like just a single spike.

figure
surfc(x,y,I)
xlim([-8 8]),ylim([-8 8])

I don't know why you would want to use such a wide range that is some 100X the actual range. Maybe a log scaling might work, let's see, don't know that I've ever tried on a 3D...

figure
surfc(x,y,I)
xlim([-8 8]),ylim([-8 8])
hAx=gca;
hAx.XScale='log'; hAx.YScale='log';

Oh yeah, ya' can't do that with negative values...but it does show the positive side

You could set the limits to something more reasonable and still see something...

figure
surfc(x,y,I)
M=0.1;
xlim(M*[-1 1]),ylim(M*[-1 1])

Warning: Negative limits ignored

Experiment with "M" -- to show anything else on the Z origin plane, you'll have to define a Z array over the larger range and then insert the finer grid in the middle...or use hold to keep the original figure and then plot the zero plane separately.

Answer 1

[dpb]

To amplify on previous to incorporate the larger range and still compute over a fine mesh in ROI. I converted the former Answer back to a Comment after realized the simple way to do it later...

The following also removes the references to the Symbolic TB--it isn't needed here if don't try to evaluate exponentials at absurd argument values and expect to get finite results.

You'll see the fine mesh in the middle with the large number of grid lines; it's smooth in the X dimension such that the grid there could be quite a lot coarser with no loss of information; Y is the culprit...

But, this builds the plane at the Z origin as a single array although again why to use such a large region is mysterious...

x0 = -0.03:0.001:0.03;
x1=[-8 -4 -1 -0.05];                % augment to cover the large range coarsely
x0=[x1 x0 sort(abs(x1))];           % add to fine range front and back...
y0 = x0;
[x,y] = meshgrid(x0) ;
z=1;
lambda=532E-9;
wo = 0.01;
k = (2.*pi)./(lambda);
Pho=0.0750;
Gamma     = ((1i.*k)./(2.*z)) + (1./(Pho.^2)) + (1./(wo.^2));
Gamma_star=((-1i.*k)./(2.*z)) + (1./(Pho.^2)) + (1./(wo.^2));
%Gamma_star = subs(Gamma, 1i, -1i);                             % don't need symbolic here
alpha = (Gamma_star) - 1./((Gamma.*Pho.^4));
beta = 1 + 1./((Pho.^2).*Gamma);
beta_s = -1 + 1./((Pho.^2).*Gamma);
C = (1./(4.*(z.^2))).*((k.^2)./(alpha.*Gamma)) ;
fx = ((1i.*k)./(2.*z)) + ((k.*k)./(4.*z.*z.*Gamma)) +  ((k.*k.*beta.*beta)./(4.*z.*z.*alpha));
fy = ((-1i.*k)./(2.*z)) + ((k.*k)./(4.*z.*z.*Gamma)) +  ((k.*k.*beta.*beta)./(4.*z.*z.*alpha));
X= C.*exp( -(fx.*x.^2) - (fy.*y.^2) ) ;
%I= abs(double((vpa(X))))./max(abs(double((vpa(X)))));
I=abs(X)./max(abs(X));
%plot(x0,I)
surfc(x,y,I)
M=0.05;                         % set limit so isn't just spike...
xlim(M*[-1 1]),ylim(M*[-1 1])
xlabel('X'),ylabel('Y')
colorbar

●

[dpb]

Illustration to make grid of appropriate resolution for each direction...

i also took some time cleaned up a bunch of the superfluous parentheses and unneeded "dot" * and / operators to reduce some of the clutter as well as eliminating the Gamma_star and beta_star variables; Gamma_star is simply the conjugate of Gamma and is only needed once in the expression for alpha so just added the transpose operator ' there in line. beta_star was never used so simply removed it. That's all basically just aesthetics, but makes the code easier to read and removes a few unneeded calculations besides along the way.

dx=0.01; dy=0.001;                  % use variables, don't bury magic numbers in code
xRnge=0.3; yRnge=0.3;               % could also make different lower/upper if chose
x0=-xRnge:dx:xRnge;                 % base x 
x1=[-8 -4 -1 -0.05];                % augment to cover the large range coarsely
x0=[x1 x0 sort(abs(x1))];           % add to fine range front and back...
y0=-yRnge:dy:yRnge;                 % ditto for y now...
y0=[x1 y0 sort(abs(x1))];           % add to fine range front and back...
[x,y] = meshgrid(x0,y0);            % use different grid each direction
z=1;
lambda=532E-9;
wo = 0.01;
k = 2*pi/lambda;
Pho=0.0750;
Gamma=1i*k/(2*z) + 1/Pho^2 + 1/wo^2;
alpha=Gamma'-1/(Gamma*Pho.^4);
beta  =  1 + 1/(Pho^2*Gamma);
C = (1/(4*z^2))*(k^2/(alpha.*Gamma));
fx = ( 1i*k)/(2*z) + (k*k)/(4*z*z*Gamma) + (k*k*beta*beta)/(4*z*z*alpha);
fy = (-1i*k)/(2*z) + (k*k)/(4*z*z*Gamma) + (k*k*beta*beta)/(4*z*z*alpha);
X= C.*exp(-(fx.*x.^2) - (fy.*y.^2)) ;
I=abs(X)./max(abs(X));
surfc(x,y,I)
Mx=0.1; My=0.05;                         % set limit so isn't just spike...
xlim(Mx*[-1 1]),ylim(My*[-1 1])
xlabel('X'),ylabel('Y')
colorbar

●

[Athira T Das]

@dpb Thank you for the clarification.

As to make the contour plot visible fully, I need to lift the surf plot. How it can be done?

Some thing like as shown in the attached figure.

[dpb]

dx=0.01; dy=0.001;                  % use variables, don't bury magic numbers in code
xRnge=0.3; yRnge=0.3;               % could also make different lower/upper if chose
x0=-xRnge:dx:xRnge;                 % base x 
x1=[-8 -4 -1 -0.05];                % augment to cover the large range coarsely
x0=[x1 x0 sort(abs(x1))];           % add to fine range front and back...
y0=-yRnge:dy:yRnge;                 % ditto for y now...
y0=[x1 y0 sort(abs(x1))];           % add to fine range front and back...
[x,y] = meshgrid(x0,y0);            % use different grid each direction
zRnge=[-0.5 1];                     % the z axis range for display
z=1;
lambda=532E-9;
wo = 0.01;
k = 2*pi/lambda;
Pho=0.0750;
Gamma=1i*k/(2*z) + 1/Pho^2 + 1/wo^2;
alpha=Gamma'-1/(Gamma*Pho.^4);
beta  =  1 + 1/(Pho^2*Gamma);
C = (1/(4*z^2))*(k^2/(alpha.*Gamma));
fx = ( 1i*k)/(2*z) + (k*k)/(4*z*z*Gamma) + (k*k*beta*beta)/(4*z*z*alpha);
fy = (-1i*k)/(2*z) + (k*k)/(4*z*z*Gamma) + (k*k*beta*beta)/(4*z*z*alpha);
X= C.*exp(-(fx.*x.^2) - (fy.*y.^2)) ;
I=abs(X)./max(abs(X));
surfc(x,y,I)
Mx=0.1; My=0.05;                    % set limit so isn't just spike...
xlim(Mx*[-1 1]),ylim(My*[-1 1])
xlabel('X'),ylabel('Y')
colorbar
zlim(zRnge)
hAx=gca;                            % get handle to axes
hAx.ZAxis.TickLabelFormat='%0.1f';  % make uniform formatting
hAx.ZTick=hAx.ZTick(hAx.ZTick>=0);  % don't show negative ticks

●

Salt to suit... :)

[Athira T Das]

@dpb Great!!!!

[Athira T Das]

@dpb I tried to plot with different I value. Attaching the code. But the still the plot is not clearly visible. Can you help what is wrong in the code?

clc;close all;clear all;
dx=0.01; dy=0.001;                  % use variables, don't bury magic numbers in code
xRnge=0.3; yRnge=0.3;               % could also make different lower/upper if chose
x0=-xRnge:dx:xRnge;                 % base x 
x1=[-8 -4 -1 -0.05];                % augment to cover the large range coarsely
x0=[x1 x0 sort(abs(x1))];           % add to fine range front and back...
y0=-yRnge:dy:yRnge;                 % ditto for y now...
y0=[x1 y0 sort(abs(x1))];           % add to fine range front and back...
[x,y] = meshgrid(x0,y0);   
zRnge=[-0.5 1];                     % the z axis range for display
z=100;
m=1;Omega=0.1;sgm=1;
lambda=532E-9;
wo = 0.01;
k = 2*pi/lambda;
ettaa=10^-3;
chiT=10^-8;
w=-2.5;
epsilon=10^-5;
Cn=(1/4.6416)*(chiT*(epsilon)^(-1/3))*(10^-8) *(10^2)
fun1=1.28.*(k.^2).*z.*((ettaa).^(-1/3)).*(Cn).*(6.78 + (47.57./(w.^2)) - (17.67./w) ) ;
Pho =fun1.^(-0.5)


 u=1
    Gamma=1i.*k./(2.*z) + 1./Pho.^2 + 1./wo.^2;
    alpha=Gamma'-1./(Gamma.*Pho.^4);
    C = (1./(4.*(lambda.^2).*(z(u).^4))).*((pi.^2)./(alpha.*Gamma)).*(1./(2.*1i.*sqrt(Gamma)).^m); 
    fx = ((Omega.^2)./(4.*Gamma)) + ((1i.*k.*x.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*x.^2)./(4.*z(u).^2.*Gamma));
    fy = ((Omega.^2)./(4.*Gamma)) + ((1i.*k.*y.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*y.^2)./(4.*z(u).^2.*Gamma));
    gx = ((Omega.^2)./(4.*Gamma)) - ((1i.*k.*x.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*x.^2)./(4.*z(u).^2.*Gamma));
    gy = ((Omega.^2)./(4.*Gamma)) - ((1i.*k.*y.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*y.^2)./(4.*z(u).^2.*Gamma));
    Bx_p = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
    Bx_n = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
    By_p = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
    By_n = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
    Kx_p = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
    Kx_n = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
    Ky_p = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
    Ky_n = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
    
    X = 0;
    for l = 0:m
      L = (((1i).*sgm).^(m-l)).*nchoosek(m,l);
      sum_j = 0;
      for j = 0:m
          J= (((-1i).*sgm).^(m-j)).*nchoosek(m,j);  
          
          sum1 = 0.0;
          for p = 0:1/2:l/2
                faktor1_P = (((-1).^p).*factorial(l))./(gamma(p+1).*factorial(l-(2.*p))).*(((2.*1i)./(sqrt(Gamma))).^(l-(2.*p))).*exp(fx).*exp(fy);
                for s1 = 0:l-2*p
                    faktor1_s1 =  nchoosek((l-(2.*p)),s1).*((1./(Pho(u).^2)).^s1).*((((1i.*k.*x)./(2.*z(u)))+(Omega./2)).^(l-(2.*p)-s1)).*(1./(((2.*1i.*sqrt(alpha))).^(j+s1)));
                    for q = 0:1/2:(m-l)/2
                        faktor1_Q = (((-1).^q).*factorial(m-l))./(gamma(q+1).*factorial(m-l-(2.*q))).*(((2.*1i)./(sqrt(Gamma))).^(m-l-(2.*q)));
                        sum_inner_1 = 0;
                        for s2 = 0:m-l-2*q
                            E1 = (exp(((Bx_p).^2)./alpha).*hermiteH(j+s1,((1i.*Bx_p)./(sqrt(alpha)))).*exp(((By_p).^2)./alpha).*hermiteH(m-j+s2,((1i.*By_p)./(sqrt(alpha))))) - (exp(((Bx_n).^2)./alpha).*hermiteH(j+s1,((1i.*Bx_n)./(sqrt(alpha)))).*exp(((By_n).^2)./alpha).*hermiteH(m-j+s2,((1i.*By_n)./(sqrt(alpha)))));
                            sum_inner_1 = sum_inner_1 + nchoosek((m-l-(2.*q)),s2).*((1./(Pho(u).^2)).^s2).*((((1i.*k.*y)./(2.*z(u)))+(Omega./2)).^(m-l-(2.*q)-s2)).*(1./(((2.*1i.*sqrt(alpha))).^(m-j+s2))).*E1;
                        end
                    end
                end
                sum1 = sum1 + faktor1_P.*faktor1_s1.*faktor1_Q.*sum_inner_1;
           end
      
           
       sum_j= sum_j + J.*(sum1);
      end
     X = X + (L.*sum_j);
    end  
    str2 = num2str(z(u));
    str1 = 'SGV_z=';
    tname1 = strcat(str1,str2,'.tif');
    
    X  = (C.*X);  
   
    I=abs(X)./max(abs(X)); 
     plot(x0,I)
    g  = figure;                                                  
    surfc(x,y,I,'FaceAlpha',1,'EdgeColor','none')
    Mx=0.1; My=0.05;                    % set limit so isn't just spike...
    xlim(Mx*[-1 1]),ylim(My*[-1 1])
    xlabel('X'),ylabel('Y')
    colorbar;
    colormap default;   
    axis('equal');axis off; 
    title(sprintf('z = %.2f m', z(u)));
    exportgraphics(g,tname1);

[dpb]

...    
surfc(x,y,I,'FaceAlpha',1,'EdgeColor','none')
Mx=0.1; My=0.05;                    % set limit so isn't just spike...
xlim(Mx*[-1 1]),ylim(My*[-1 1])

If you're just continuing to use pieces of the prior Answer that were determined empirically to be useful in the prior case with new data, it's not at all surprising they may not work well for some other set of data.

I tried to run the above code; it took over 4 minutes and time out here and I've got too much else going on in the workspace to be able to mung on this on local machine at the moment, sorry.

You just need to again start at the basics of ensurig what you've calculated is what you intended and correct and then figure out the way in which to display the plot effectively utilizing the various orienations and scaling tools to do so.

I'd note the original is still full of a lot of NaN values that are simply hidden by the final plot; it's quite possible something like that is also happening in your new case. It may not hurt to do as I did to get started before; cut down the size to something much more manageable to begin with and then debug from there...having no idea what it is you're trying to do makes doing anything about the actual code itself almost impossible besides it again being so replete with redundant parentheses and repeated terms that it is almost impossible to read and get any idea of what is actually trying to be computed.

[dpb]

>> athira
'hermiteH' requires Symbolic Math Toolbox.
Error in athira (line 55)
              E1 = (exp(((Bx_p).^2)./alpha).*hermiteH(j+s1,((1i.*Bx_p)./(sqrt(alpha)))).*exp(((By_p).^2)./alpha).*hermiteH(m-j+s2,((1i.*By_p)./(sqrt(alpha))))) - (exp(((Bx_n).^2)./alpha).*hermiteH(j+s1,((1i.*Bx_n)./(sqrt(alpha)))).*exp(((By_n).^2)./alpha).*hermiteH(m-j+s2,((1i.*By_n)./(sqrt(alpha))))); 
>>

I don['t have the Symbolic TB, either.

That line is over 300(!!!) characters in length; trying to debug anything like that is essentially impossible...

hermiteH

[dpb]

I converted to a function so wouldn't pollute/destroy my working environment and removed the reference to the Symbolic TB hermiteH function with a <submission from the FEX> instead...

Reverting back to the same case as before, simply making the ranges small enough to see what is there; it doesn't look bad; rotating it has some pretty interesting views...I've absolutely no klew what these must be representing, but you just need to look more closely at what you're doing in code with the viewpoints and scaling -- experiment at the command line first to see what is going on -- things don't look bad until you added at the very end I commented back out...

%function athira
  dx=0.01; dy=0.001;                  % use variables, don't bury magic numbers in code
  xRnge=0.3; yRnge=0.3;               % could also make different lower/upper if chose
  x0=-xRnge:dx:xRnge;                 % base x 
  x1=[-8 -4 -1 -0.05];                % augment to cover the large range coarsely
  x0=[x1 x0 sort(abs(x1))];           % add to fine range front and back...
  y0=-yRnge:dy:yRnge;                 % ditto for y now...
  y0=[x1 y0 sort(abs(x1))];           % add to fine range front and back...
  [x,y] = meshgrid(x0,y0);   
  zRnge=[-0.5 1];                     % the z axis range for display
  z=100;
  m=1;Omega=0.1;sgm=1;
  lambda=532E-9;
  wo = 0.01;
  k = 2*pi/lambda;
  ettaa=10^-3;
  chiT=10^-8;
  w=-2.5;
  epsilon=10^-5;
  Cn=(1/4.6416)*(chiT*(epsilon)^(-1/3))*(10^-8) *(10^2);
  fun1=1.28.*(k.^2).*z.*((ettaa).^(-1/3)).*(Cn).*(6.78 + (47.57./(w.^2)) - (17.67./w) ) ;
  Pho =fun1.^(-0.5);
  u=1;
    Gamma=1i.*k./(2.*z) + 1./Pho.^2 + 1./wo.^2;
    alpha=Gamma'-1./(Gamma.*Pho.^4);
    C = (1./(4.*(lambda.^2).*(z(u).^4))).*((pi.^2)./(alpha.*Gamma)).*(1./(2.*1i.*sqrt(Gamma)).^m); 
    fx = ((Omega.^2)./(4.*Gamma)) + ((1i.*k.*x.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*x.^2)./(4.*z(u).^2.*Gamma));
    fy = ((Omega.^2)./(4.*Gamma)) + ((1i.*k.*y.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*y.^2)./(4.*z(u).^2.*Gamma));
    gx = ((Omega.^2)./(4.*Gamma)) - ((1i.*k.*x.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*x.^2)./(4.*z(u).^2.*Gamma));
    gy = ((Omega.^2)./(4.*Gamma)) - ((1i.*k.*y.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*y.^2)./(4.*z(u).^2.*Gamma));
    Bx_p = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
    Bx_n = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
    By_p = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
    By_n = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
    Kx_p = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
    Kx_n = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
    Ky_p = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
    Ky_n = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
    
    X = 0;
    for l = 0:m
      L = (((1i).*sgm).^(m-l)).*nchoosek(m,l);
      sum_j = 0;
      for j = 0:m
        J= (((-1i).*sgm).^(m-j)).*nchoosek(m,j);  
        sum1 = 0.0;
        for p = 0:1/2:l/2
          faktor1_P = (((-1).^p).*factorial(l))./(gamma(p+1).*factorial(l-(2.*p))).*(((2.*1i)./(sqrt(Gamma))).^(l-(2.*p))).*exp(fx).*exp(fy);
          for s1 = 0:l-2*p
            faktor1_s1 =  nchoosek((l-(2.*p)),s1).*((1./(Pho(u).^2)).^s1).*((((1i.*k.*x)./(2.*z(u)))+(Omega./2)).^(l-(2.*p)-s1)).*(1./(((2.*1i.*sqrt(alpha))).^(j+s1)));
            for q = 0:1/2:(m-l)/2
              faktor1_Q = (((-1).^q).*factorial(m-l))./(gamma(q+1).*factorial(m-l-(2.*q))).*(((2.*1i)./(sqrt(Gamma))).^(m-l-(2.*q)));
              sum_inner_1 = 0;
              for s2 = 0:m-l-2*q
                E1 = (exp(((Bx_p).^2)./alpha).*hermiteH(j+s1,(1i.*Bx_p)/sqrt(alpha)).*exp(((By_p).^2)./alpha).*hermiteH(m-j+s2,((1i.*By_p)./(sqrt(alpha))))) - (exp(((Bx_n).^2)./alpha).*hermiteH(j+s1,((1i.*Bx_n)./(sqrt(alpha)))).*exp(((By_n).^2)./alpha).*hermiteH(m-j+s2,((1i.*By_n)./(sqrt(alpha)))));
                sum_inner_1 = sum_inner_1 + nchoosek((m-l-(2.*q)),s2).*((1./(Pho(u).^2)).^s2).*((((1i.*k.*y)./(2.*z(u)))+(Omega./2)).^(m-l-(2.*q)-s2)).*(1./(((2.*1i.*sqrt(alpha))).^(m-j+s2))).*E1;
              end
            end
          end
          sum1 = sum1 + faktor1_P.*faktor1_s1.*faktor1_Q.*sum_inner_1;
        end
        sum_j= sum_j + J.*(sum1);
       end
       X = X + (L.*sum_j);
    end  
    str2 = num2str(z(u));
    str1 = 'SGV_z=';
    tname1 = strcat(str1,str2,'.tif');
    X  = (C.*X);  
    I=abs(X)./max(abs(X)); 
  %  plot(x0,I)
  %  hF  = figure;                                                  
    surfc(x,y,I,'FaceAlpha',1,'EdgeColor','none')
    Mx=0.1; My=0.05;                    % set limit so isn't just spike...
    xlim(Mx*[-1 1]),ylim(My*[-1 1])
    zlim(zRnge)
    xlabel('X'),ylabel('Y')
    hAx=gca;
    set([hAx.XAxis;hAx.YAxis;hAx.ZAxis],{'TickLabelFormat'},{'%0.2f';'%0.2f';'%0.1f'})
    ztk=zticks; zticks(ztk(ztk>=0))
    hAx.CameraPosition=[-1.5 0.15 5];

  %  colorbar;
  %  colormap default;   
  %  axis('equal');axis off; 
  %  title(sprintf('z = %.2f m', z(u)));
  %  exportgraphics(hF,tname1);
%end

function res=hermiteh(n,x)
% Useage:
%  hermiteh(n,x)
%     generate the nth hermite polynomial of x

  m=0:floor(n/2);
  q2=width(m);
  s=ndims(x);
  [g1,g2]=size(x);
  X=repmat(x,[ones(1,s), q2]);
  m=permute(repmat(m,[g1,1,g2]),[1,3,2]);
  res=factorial(n)*sum((-1).^m./(factorial(m).*factorial(n-2*m)).*(2*X).^(n-2*m),3);
end

function h=hermiteH(n,x)
% replace calls to Symbolic TB with FEX submission...
  h=hermiteh(n,x);
end