How do I plot a solution from a PDE within a unit square?

I have solved a nonlinear PDE but now need to plot the solution as a line plot rather than as a surface plot. The solution matrix is given as (phi0) and the unit square is of x and y dimension. How do I obtain a line plot of the solution within this unit square ? see attached for something similar of what i intend doing.

clear;

clc

% initialise number of nodes

L=1;

H=1;

M=10; % number of x nodes

N=10;

dx=L/(M-1);

dy=H/(N-1);

x=linspace(0,L,M); % discreet x values

y=linspace(0,H,N);

% initialise Guess

phi0=ones(M,N)*0.5

%Boundary conditions phi0(1,:)=0; % bottom boundary

phi0(M,:)=1;% right boundary

phi0(:,1)=0; %left boundary

phi0(:,N)=1; % top boundary

% initialise Guess

phi_delta=ones(M,N)*0.5 %zeros(M,N);

% initialise Boundary conditions

phi_delta0(1,:)=0; % bottom boundary

phi_delta0(M,:)=1;% right boundary

phi_delta0(:,1)=0; %left boundary

phi_delta0(:,N)=1; % top boundary

Tol=10^-6;

Tol1=10^-6;

flag1=2

while flag1==2

    flag1=0;

    flag = 1;

    while flag==1

        flag =0;

        phi_delta=phi_delta0

        for i=2:M-1

            for j=2:N-1

                f(i,j)=((phi0(i+1,j)^2-phi0(i-1,j)^2)/2*dx)+((phi0(i,j+1)^2-phi0(i,j-1)^2)/2*dy)-((phi0(i+1,j)-2*phi0(i,j)+phi0(i-1,j)/dx^2)+...
                    (phi0(i,j+1)-2*phi0(i,j)+phi0(i,j+1)/dy^2))-1/2*((phi0(i+1,j)^2-2*phi0(i,j)^2+phi0(i-1,j)^2/dx^2)...
                    +(phi0(i,j+1)^2-2*phi0(i,j)^2+phi0(i,j-1)^2/dy^2));

                aE=(phi0(i+1,j)/dx)-((1+phi0(i+1,j))/dx^2);

                aW=(-phi0(i-1,j)/dx)-((1+phi0(i-1,j))/dx^2);

                aN=(phi0(i,j+1)/dy)-((1+phi0(i,j+1))/dy^2);

                aS=(-phi0(i,j-1)/dy)-(1+phi0(i,j-1))/dy^2;

                a0=(2*(1+phi0(i,j))/dx^2)+(2*(1+phi0(i,j))/dy^2); 

                phi_delta0(i,j)=(-f(i,j)-aE*phi_delta(i+1,j)-aW*phi_delta(i-1,j)-aN*phi_delta(i,j+1)-aS*phi_delta(i,j-1))/a0

                error = abs((phi_delta0(i,j)-phi_delta(i,j))/phi_delta(i,j));

                if error > Tol

                    flag=1

                end
            end

        end

    end

    phi0=phi0+0.05*(phi_delta0)

    if error > Tol1

        flag1=2

    end

end

%Boundary conditions

phi0(1,:)=0; % bottom boundary

phi0(M,:)=1;% right boundary

phi0(:,1)=0; %left boundary

phi0(:,N)=1; % top boundary

Answer 1

[John D'Errico]

help contour

[JohnF]

Thanks john, i tried before using contour(x,y,phi0) but what i get is multiple lines within the square instead of a single line representing the solution. I know the pic i showed had multiple lines but that's because it was for various convergence criteria.