Lorenz Equation using Newton's Method

I am doing my project on writing Matlab code for the Lorenz equation using Newton's Method. My task was to write a code by using while loop so that the roots converge. I have posted my code below, where I couldn't able get the convergence.

r=28; sigma=10; beta=8/3; 
x1=0; y1=0; z1=0;
x2=sqrt(beta*(r-1)); y2=sqrt(beta*(r-1)); z2=r-1;
x3=-sqrt(beta*(r-1)); y3=-sqrt(beta*(r-1)); z3=r-1;
nx=500; nz=500;
xmin=-40; xmax=40; zmin=-40; zmax=40;
x_grid=linspace(xmin,xmax,nx); z_grid=linspace(zmin,zmax,nz);
[X,Z]=meshgrid(x_grid,z_grid);
RelTol=1.e-06; AbsTol=1.e-09;

            for i=1:3
        if i==1 , x=x1; y=y1; z=z1; end
        if i==2 , x=x2; y=y2; z=z2; end
        if i==3 , x=x3; y=y3; z=z3; end
        error=Inf;
        for j=1:nx
              for k=1:nz
                y0=3*sqrt(2);
                while error<=max(RelTol*max(abs([x,y,z])),AbsTol)
                    J  =   [-sigma, sigma,0;r-z_grid(k),-1,-x_grid(j);y0,x_grid(j),-beta];
                   rhs =  -[(sigma*(y0-x_grid(j)));(x_grid(j)*(r-z_grid(k))-y0);((x_grid(j)*y0)-(beta*z_grid(k)))];
                   delta_xyz= J\rhs;
                    x_grid(j) = x_grid(j) + delta_xyz(1);
                    y0 = y0+delta_xyz(2);
                    z_grid(k) = z_grid(k) + delta_xyz(3);
                    error=max(abs(delta_xyz));
                end 
            X(j,k)=x_grid(j);
            Z(k,j)=z_grid(k);
              end
        end
end