clc
V_t=26e-3;
I_s=1e-6;
I_d=0.1
V_d=V_t*log(1 + I_d/I_s)
G_eq=I_s/V_t*exp(V_d/V_t)
I_eq=I_d- G_eq*V_d
I_dold = 1;
R = [1 2 1];
V = [1 0 1];
d = [2 0];
abs(I_dold - I_eq)

c = 1;

n = max([max(R(:, 1: 2)) max(V(:, 1:2)) max(D)]);
m = size(V, 1);
while(abs(I_dold-I_d)>1e-1)
A = zeros(n + m);
G = zeros(n);
B = zeros(n, m);
C = zeros(m, n);
D = zeros(m, m);

for i = 1 : size(R, 1)
   n1 = R(i, 1)
   n2 = R(i, 2)
   g = 1/R(i, 3)
   if(n1 == 0)
       G(n2, n2) = G(n2, n2) + g
   elseif (n2 == 0)
       G(n1, n1) = G(n1, n1) + g
   else
       G(n1, n2) = G(n1, n2) - g
       G(n2, n1) = G(n2, n1) - g
       G(n2, n2) = G(n2, n2) + g
       G(n1, n1) = G(n1, n1) + g
   end
end
G
for i = 1 : size(d, 1)
   np = d(i, 1);
   nn = d(i, 2);
   g  = G_eq;
   if(np == 0)
       G(nn, nn) = G(nn, nn) + g
   elseif (nn == 0)
       G(np, np) = G(np, np) + g
   else
       G(np, nn) = G(np, nn) - g
       G(nn, np) = G(nn, np) - g
       G(nn, nn) = G(nn, nn) + g
       G(np, np) = G(np, np) + g
   end

end
G


for i = 1 : size(V, 1)
   nplus = V(i, 1);
   nminus = V(i, 2);
   if (nplus == 0)
       B(nminus, i) = B(nminus, i) - 1;
   elseif (nminus == 0)
       B(nplus, i) = B(nplus, i) + 1;
   else
       B(nminus, i) = B(nminus, i) - 1;
       B(nplus, i) = B(nplus, i) + 1;
   end
end
B

C = B';

A = [G B; C D]

z = zeros(n + m, 1);
i = zeros(n, 1);
e = zeros(m, 1);

for j = 1 : size(V, 1)
   e(j) = V(j, 3);
end

for k = 1 : size(i, 1)
   for j = 1 : size(d, 1)
      if d(j, 1) == k
          i(k) = i(k) - I_eq;
      elseif d(j, 2) == k
          i(k) = i(k) + I_eq;
      end
   end
end

z = [i; e]

X = inv(A)*z

I_dold = I_eq
I_eq = -X(end);
V_d=V_t*log(1 + I_eq/I_s)
G_eq=I_s/V_t*exp(V_d/V_t)
end