Vd1 : 0.216 + 0.001 $ Id1 : 0.1e-3 $  /* 実測データその１, Vd を少し大きめに補正 */
Vd2 : 0.341 - 0.02  $ Id2 : 10e-3  $  /* 実測データその２, 10mA * 2Ω = 0.02 */

Is : 1e-7 $ N : 1 $ VT : 0.026 $

for i : 0 while i < 100 do (
  /* Is を漸近させる */
  tId1 : Is*(exp(Vd1/(N*VT)) - 1),
  tId2 : Is*(exp(Vd2/(N*VT)) - 1),
  u1 : tId1/Id1,
  u2 : tId2/Id2,
  if abs(1 - u1) < 1e-6 and abs(1 - u2) < 1e-6 then return(i),  /* 終了？ */
  Is : Is / (u1*u2)**0.4,

  /* N を漸近させる */
  tId1 : Is*(exp(Vd1/(N*VT)) - 1),
  tId2 : Is*(exp(Vd2/(N*VT)) - 1),
  u1 : tId1/Id1,
  u2 : tId2/Id2,
  N : N * (u2/u1)**0.25  /* 18@0.25 */
);
display(Is, N, tId1, tId2) $