Relaxation Example 3 MATLAB

%Clear command window, saved variables, and any active figures.
clc
clear
clf
%Parameters
phi1=200;
phi2=-phi1;
iterations=5000;
%Outer length, width, and grid spacing in centimeters
len=300;
sep=5;
K=len/sep+1;

%Initialize the array.
phi=zeros(K,K);

%Position on grid of left end of outlet and its length
X_out = 200/sep+1;
Y_out = 1;
out_len = 100/sep-1;

%Position on grid of upper end of inlet and its length
X_in = 1;
Y_in = 200/sep+1;
in_len=100/sep-1;

%Impose boundary conditions
for k = X_out:(X_out+out_len)
    phi(k,Y_out) = phi1;
end
for k = Y_in:(Y_in+in_len)
    phi(X_in,k) = phi2;
end
%Now carry out the relaxation.
for iter=1:iterations
    for i=1:K
        for j=1:K
            if i==K || j==K || i==1 || j==1
                phi(i,j) = phi(i,j);
            elseif i==(X_out-1) && j<Y_in
                phi(i,j) = phi(i,j);
            elseif j==(Y_in-1) && i<X_out
                phi(i,j) = phi(i,j);
            else
                phi(i,j)=0.25*(phi(i+1,j)+phi(i-1,j)+phi(i,j+1)+phi(i,j-1));
            end
        end
    end
end

%We only want to graph the flow in the pipe so we'll set the grid elements
%outside the pipe to NaN
for i=1:(X_out-2);
    for j=1:(Y_in-2);
        phi(i,j) = NaN;
    end
end
phi=flip(phi');
[x y] = meshgrid(0:sep:len);
[gx gy] = gradient(phi);
startx = 200:10:300;
starty = 100:-10:0;
%Uncomment this if you want a contour and/or streamline plot
% hold on
% contourf(x,y,phi);
% streamline(x,y,gx,gy,startx,starty)
% streamline(x,y,-gx,-gy,startx, starty)
% hold off
mesh(x,y,phi)
zlabel('Potential in m^2/s')
xlabel('x')
ylabel('y')