Why contourf and fcontour gives different images for the same matlab code?

I have code for fcountour as;

syms x y  
m =1;
sgm =1;
wo = 0.01;
lambda=532 *10^-9;
k=2*pi/lambda;
Omega = .2; 
 
z=30;
Pho=0.003;


for u=1:length(z)
    I=0;
    Gamma = ((1i.*k)./(2.*z(u))) + (1./(Pho(u).^2)) + (1./(wo.^2));
    Gamma_star = subs(Gamma, 1i, -1i);
    alpha = (Gamma_star) - 1./((Gamma.*Pho(u).^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  
 
    X  = (C.*X);  
    I = I + abs(vpa(X)) ;
    figure;
    fcontour(real(I),  [-1 1  -1 1]*1E-1, 'Fill','on', 'MeshDensity',150);
    colormap('hot');  
    axis('equal');colorbar('vert');
     
end

Then i have contourf plot as;

figure;

m =1;
sgm =1;
wo = 0.1;
lambda=532 *10^-9;
k=2*pi/lambda;
Omega = .2; 
 
z=30

z = 30

Pho=0.003

Pho = 0.0030

[x, y] = meshgrid(-1:0.06:1);


for u=1:length(z)
    I=0;
    Gamma = ((1i.*k)./(2.*z(u))) + (1./(Pho(u).^2)) + (1./(wo.^2));
    Gamma_star = subs(Gamma, 1i, -1i);
    alpha = (Gamma_star) - 1./((Gamma.*Pho(u).^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  
    
    X  = (C.*X);  
    I = I + abs(vpa(X)) ;
    contourf(x,y,real(abs(I))/max(real(abs(I(:)))))
    colormap('hot');
    colorbar('vert');
    axis('equal') 
end

According to me fcountour plot is correct, then what is wrong with countourf plot? Can anyone help me to fix the error in contourf plot?

[VBBV]

What I got is plot with two extra lobes and the figure is not complete. I have done the contourf plot is to normalize my data. Why when I normalize the data, looks different?

The figures produced by both contourf & fcontour functions are same. If you look at fcontour the range limits

fcontour(real(I),  [-1 1  -1 1]*1E-1, 'Fill','on', 'MeshDensity',150);

are [ -1 1 -1 1 ] multiplied by (1/10) i.e. 0.1 and fcontour doesnt normalize data when it plots. However, when you normalize the data ( as in case of contourf) it divides with maximum value from the computed vector. Also, the absoulute value was taken, this also means, the resultant value (which includes complex or imaginary values if any) is considered instead of actual value,

Answer 1

[VBBV]

If you change the meshgrid limits in the below line

[x, y] = meshgrid(-1/10:0.005:1/10);

and de-normalize the I variable in the contourf function , you will get similar figure as fcontour

contourf(x,y,real(I))

[Athira T Das]

m =1;
sgm =1;
wo = 0.1;
lambda=532 *10^-9;
k=2*pi/lambda;
Omega = .2; 
 
z=30

z = 30

Pho=0.003

Pho = 0.0030

[x, y] = meshgrid(-1/10:0.005:1/10);


for u=1:length(z)
    I=0;
    Gamma = ((1i.*k)./(2.*z(u))) + (1./(Pho(u).^2)) + (1./(wo.^2));
    Gamma_star = subs(Gamma, 1i, -1i);
    alpha = (Gamma_star) - 1./((Gamma.*Pho(u).^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  
    
    X  = (C.*X);  
    I = I + abs(vpa(X)) ;
     contourf(x,y,real(abs(I))/max(real(abs(I(:)))), 'EdgeColor','none')
    colormap('hot');
    colorbar('vert');
    axis('equal') 
end

What I got is plot with two extra lobes and the figure is not complete. I have done the contourf plot is to normalize my data. Why when I normalize the data, looks different?

[VBBV]

m =1;
sgm =1;
% this parameter is also something different in both codes
wo = 0.01;
lambda=532 *10^-9;
k=2*pi/lambda;
Omega = .2; 
 
z=30

z = 30

Pho=0.003

Pho = 0.0030

[x, y] = meshgrid(-1/10:0.005:1/10);


for u=1:length(z)
    I=0;
    Gamma = ((1i.*k)./(2.*z(u))) + (1./(Pho(u).^2)) + (1./(wo.^2));
    Gamma_star = subs(Gamma, 1i, -1i);
    alpha = (Gamma_star) - 1./((Gamma.*Pho(u).^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  
    
    X  = (C.*X);  
    I = I + abs(vpa(X)) ;
    contourf(x,y,real((I)), 'EdgeColor','none')
    colormap('hot');
    colorbar('vert');
    axis('equal') 
end

In addition to x & y range limits in meshgrid and denormailized I variable, the value of parameter wo is used differently in contourf & fcontour functions , See the below line in both codes

% this value is 0.01 in fcontour and 0.1 in contourf function codes, why ?
wo = 0.01;