Why is the contour plot flipped i.e. the value that should be at the bottom right is plotted on the top left?

Nx = input('Number of nodes in x-direction');
Ny = input('Number of nodes in y-direction');
dx = 1/(Nx - 1);
dx2 = dx*dx;
dy = 1/(Ny - 1);
dy2 = dy*dy;
x = 0:dx:1;
y = 0:dy:1;

for i = 1:Nx
    for j = 1:Ny
      phi_an(i,j) = 500*exp(-50*((1-x(i))^2 + y(j)^2)) + 100*x(i)*(1-y(j));
    end
end


Sphi = zeros(Nx,Ny);
for i = 1:Nx
    for j = 1:Ny
      Sphi(i,j) = 50000*exp(-50*((1-x(i))^2 + y(j)^2))*(100*((1-x(i))^2 + y(j)^2)-2);
    end
end

% Defining solution vector with intial guess
phi = zeros(Nx,Ny);
phi(1,1:Ny) = 500*exp(-50*(1+y.^2)); % Left BC
phi(Nx,1:Ny) = 100*(1-y) + 500*exp(-50*y.^2); % Right BC
phi(1:Nx,1) = 100*x + 500*exp(-50*((1 - x).^2)); % Bottom BC
phi(1:Nx,Ny) = 500*exp(-50*((1-x).^2 + 1)); % Top BC

% Setting up the link coefficients

ae = 1/dx2;
aw = 1/dx2;
an = 1/dy2;
as = 1/dy2;
ao = -(2/dx2 + 2/dy2);

toc


%% Jacobi Method

tic
% Matrix to store new values

phin = zeros(Nx-2,Ny-2);

%Setting up the iterations
mat = zeros(1,1000000);
R = zeros(Nx-2,Ny-2);

for n = 1:1000000
    R2 = 0;
    for i = 2:Nx - 1
        for j = 2:Ny - 1
        phin(i-1,j-1) = (Sphi(i,j) - ae*phi(i-1,j) - aw*phi(i+1,j) - an*phi(i,j+1) - as*phi(i,j-1))/ao; 
        end
    end
    
    for r = 2:Nx-1
        for s = 2:Ny-1
        R(r-1,s-1) = Sphi(r,s) - ae*phi(r-1,s) - aw*phi(r+1,s) - an*phi(r,s+1) - as*phi(r,s-1) - ao*phi(r,s);
        R2 = R2 + R(r-1,s-1)*R(r-1,s-1);
        end
     end
    R2 = sqrt(R2);
    mat(n) = R2;
    if R2<10^-6
        break
    end
    phi(2:Nx-1,2:Ny-1) = phin(1:Nx-2,1:Ny-2);
    
    
end

[X,Y] = meshgrid(x,y);
contour(X,Y,phi);

The code to the above system of equations is the one above the image.

[Adam Danz]

Either use ndgrid instead of meshgrid or transpose phi or better yet, use the vectorization in KSSV's answer.

[X,Y] = ndgrid(x,y); % use ndgrid
contour(X,Y,phi);

or

[X,Y] = meshgrid(x,y);
contour(X,Y,phi.');  %transpose

Answer 1

[KSSV]

YOu need not to use loops....it is very slow. Check the below code:

Nx = 100 ; Ny = 100 ; 
x = linspace(0,1,Nx) ; 
y = linspace(0,1,Ny) ; 

[X,Y] = meshgrid(x,y) ; 
phi = 500*exp(-50*((1-X).^2+Y.^2))+100*X.*(1-Y) ; 
contour(X,Y,phi)

[Adarsh Choudhury]

Thanks It worked

[KSSV]

Thanks is accepting the answer.....