Hi,
You are trying to plot the series of Xr1, Xr2,  Xr3 without really computing to their values. Because your error tolerances are very tight. In your code, your have one small typo: X = NewtonRaphson1(X0) ;  has to be:  X = NewtonRaphson(X0)
If you wish to plot the evolutions of your computed roots, then here is a small edited part of your code:
for i = 1:length(x)
    for j = 1:length(y)
        X0 = [x(i);y(j)] ;
        % Solve the system of Equations using Newton's Method
        X(:,i) = NewtonRaphson(X0) ;  % To collect all evolved/computed roots
        % Locating the initial conditions according to error
        if norm(X-r1)<1e-8
            Xr1 = [X0 Xr1]  ;
        elseif norm(X-r2)<1e-8
            Xr2 = [X0 Xr2] ;
        else
            Xr3 = [X0 Xr3];          
        end        
    end
end
toc
warning('on') % Remove the warning off constraint
% Initialize figure
figure
set(gcf,'color','w')
plot(X(1, :), X(2,:), 'b-'), grid on; shg 
title('Basin of attraction for f(x,y) = x^2+y^2=9 and (x-2)^2+y^2 = 9 ')
Good luck.


