simax=16;
romax=6;
K = simax*romax;
V = 200;

nbc = simax*(romax-2);
[X Y] = meshgrid(1:simax,1:romax);
mdl = (Y-1)*simax + X;
mdl(2:end-1,1:end) = 0;

bc = zeros(nbc,2);
bc(:,1) = mdl(mdl~=0);

A = eye(K)*0;

B = zeros(K,1);

%Finite Difference Method contsnts in Cylindrical
delro = 0.25;
delsi = pi/8;
ronot = 0;

n=1;
for k = 1:K
    if k <= n*simax
        ll = n*simax-simax+2;   %lower limit
        ul = n*simax-1;         %upper limit
        lb = ll-1;              %lower boundary
        ub = n*simax;           %upper boundary
        ronot = n*delro;
        %inner node ~closer to origin
        Ao = (1-delro/2/ronot)/delro^2;
        %outer node ~towards edge
        Bo = (1+delro/2/ronot)/delro^2;
        %counterclockwise node and clockwise node
        Co = (delro/ronot/delsi)^2/delro^2;
        %original node ~the center node
        Do = -(2+2*(delro/ronot/delsi)^2)/delro^2;
        for k = ll:ul
        A(k,k) = Do;
        if k-simax > 0
            A(k,k-simax)=Ao;
        end
        if k+simax < K
            A(k,k+simax)=Bo;
        end
        if k+1 < K
            A(k,k+1) = Co;
        end
        if k-1 < 0
            A(k,k-1)=Co;
        end
        end
        %Boundary
        for k = ub
        Ao = 1-delro/2/ronot;
        Bo = 1+delro/2/ronot;
        Co = (delro/ronot/delsi)^2;
        Do = -(2+2*(delro/ronot/delsi)^2);
        A(k,k)=Do;
            if k-simax > 0
                A(k,k-simax)=Ao;
            end
            if k+simax < K
                A(k,k+simax)=Bo;
            end
            if k-simax+1 > 0
                A(k,k-simax+1) = Co;
            end
            if k-1 > 0
                A(k,k-1) = Co;
            end
        end
        for k = lb
        Ao = 1-delro/2/ronot;
        Bo = 1+delro/2/ronot;
        Co = (delro/ronot/delsi)^2;
        Do = -(2+2*(delro/ronot/delsi)^2);
        A(k,k)=Do;
            if k-simax > 0
                A(k,k-simax) = Ao;
            end
            if k+simax < K
                A(k,k+simax) = Bo;
            end
            if k+1 > 0
                A(k,k+1) = Co;
            end
            if k+simax-1 < K
                A(k,k+simax-1) = Co;
            end
        end 
    else
        n = n+1;
    end
end

%Boundary Conditions
for k = 1:simax
    A(k,k) = 0;
end
for k = 33:48
    A(k,k) = 0;
end
%B-Matrix
B(17:32) = -V;
x = A\B;
theta(1:96,1) = 0;
rho(1:96,1) = 0;
for g = 1:6
    up = g*simax;
    low = g*simax+1-simax;
    Ro = g*delro;
    theta(low:up)=1:16;
    rho(low:up) = Ro;
end
theta = theta*pi/8;
%x(:,3)=x;
%x(:,1)=theta;
%x(:,2)=rho;
[X Y Z] = pol2cart(theta,rho,x);
Phi = reshape(x,simax,romax);
%figure;image(Phi'); axis image; set(gcf,'color','w');
figure;contour(Phi',romax*simax); grid on; set(gcf,'color','w');