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

nbc = simax*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);

B = zeros(K,1);

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

%inner node ~closer to origin
Ao = 1-delro/2/ronot;

%outer node ~towards edge
Bo = 1+delro/2/ronot;

%counterclockwise node and clockwise node
Co = (delro/ronot/delsi)^2;

%original node ~the center node
Do = -(2+2*(delro/ronot/delsi)^2);

for k =1:K
    for n = 1:romax
        ronot=ronot+0.25;
        ul = n*simax-1;
        ll = n*simax-simax+2;
        lb = ll-1;
        ub = n*simax;
        %Regular Nodes
        for k = ll:ul
        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 < K
                A(k,k+1) = Co;
            end
            if k-1 > 0
                A(k,k-1) = Co;
            end
        end
        %Lower Boundary Nodes
        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 < K
                A(k,k+1) = Co;
            end
            if k+simax-1 < K
                A(k,k+simax-1)=Co;
            end
        end
        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
        %Boundary Conditions
        for k = 1:simax
            A(k,k+simax) = 0;
        end
        for k = 33:48
            A(k,k-simax) = 0;
        end
    end
end

%B-Matrix
B(17:32) = -V*delro^2;
x = A\B