Relaxation Example 2 MATLAB

clc
clear
%Parameters
phi1 = 10000;
phi2 = 0;
initphi = 500;
iterations = 10000;

%Length of energized plate in meters.
len=100;

%The grid spacing must be less than or equal to the thickness of the
%'bars' of the Faraday cage.
sep=1;
K=len/sep+1;

%Initialize the array. We need to pad the boundary with zeros for reasons
%explained in the comment above the relaxation stage.
phi=zeros(K+2,K+2);
for i=2:K+1
    for j=2:K+1
        phi(i,j) = initphi;
    end
end
%Separation of cage bars
bar_sep = 10;

%Impose boundary conditions
for k = 2:K+1
    phi(k,2) = phi1;
end
%At 50cm to the right, place evenly spaced bars separated by bar_sep
%Comment out this part to remove the cage.
for k = 1:bar_sep/sep:K
    phi(k, 50/sep) = 0;
end

%See the explanation in the code for Example 1. In this example, we also
%have to contend with grid points that are at the edge of the array and
%where no boundary condition is given. For these points, we average over the
%three nearest points. To avoid a mess of conditionals, we accomplish this
%by putting an extra 'layer' of zeros on the outside of the part of the
%grid we're interested in, adding up all neighboring points, but then dividing by 3 instead of 4.
%While this does mean a slight loss of accuracy at those points, the actual
%physics that we're interested in is taking place in the interior of the grid.
for iter = 1:iterations
%Uncomment this if you want to be able to track the program's progress.
%     if mod(iter, 1000) == 0
%         iter
%     end
    for i = 2:K+1
        for j = 2:K+1
            switch phi(i,j)
                case phi1
                    phi(i,j) = phi(i,j);
                case phi2
                    phi(i,j) = phi(i,j);
                otherwise
                    if i==2 || i==K+1 || j==2 || j==K+1
                        phi(i,j) = 0.33*(phi(i+1,j)+phi(i-1,j)+phi(i,j+1)+phi(i,j-1));
                    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
end
%Now we get rid of the layer of zeros to make plotting easier.
phi(:,1) = [];
phi(:,K+1) = [];
phi(1,:) = [];
phi(K+1,:) = [];
[x,y]=meshgrid(0:sep:len);
%contourf(x,y,phi);
meshc(x,y,phi);