Problem in numerical integration

Hi every one.

I want to by means of "Trapz" do numerical integration but when I plot the result data, its shape is not what I want. in my code, u1 is my function that I wanted to do integration of it. I do double integration, first is based on y1 and second is base on x1. If u1 is plotted, it is not same of u2 shape. I couldn't detect codes defect. so, who can help me?

 close all
 clear all
 clc

 N=100;
 x1=linspace(-1,1,N);
 y1=x1;
 z=0.42;
 lambda=1030e-9;
 k=2*pi/lambda;
 w0=0.1;
 zr=pi*w0.^2/lambda;
 R=z*(1+(zr/z).^2);
 w=w0*sqrt(1+(z/zr)^2);

 [x,y]=meshgrid(linspace(-1,1,N));
 [x2,y2]=meshgrid(linspace(-1,1,N));

 r=sqrt(x.^2+y.^2);
 teta=atan(y./x);

 r1=sqrt(x2.^2+y2.^2);
 teta1=atan(y2./x2);

 u1=ones(N,N);
 u=u1;
 u_2=u;
 u2=u;
 B=u;
 for i=1:N
     for j=1:N
         u(i,j)=(r(i,j)/w).*exp(-(r(i,j)/w).^2);
         B(i,j)=besselj(1,-k.*r(i,j).*r1(i,j)./z);
         u1(i,j)=r(i,j).*u(i,j).*exp(-1i*k*(r(i,j)./w).^2./(2*z)).*B(i,j);
      end
  end 

for i=1:N
    for j=1:N
      u_2(:,j)=trapz(y1,u1(:,j));
   end
   u2(i,:)=cos(teta1(i,j)).*(-k/z).*exp(-1i*k.*(r1(i,:)./w).^2./(2*z)).*trapz(x1,u_2(i,:));%
end