clc;
clear all;
close all;


global R1 R2 R3 R4 V0

V0 = 9.0;
R1 = 4700.0;
R2 = 4700.0;
R3 = 470.0;
R4 = 4700.0;

global D_Vt D_N D_Is D_Rs


#.model 1N4148 D(Is=2.52n Rs=.568 N=1.752 Cjo=4p M=.4 tt=20n Iave=200m Vpk=75 mfg=OnSemi type=silicon)

D_Is=2.52e-9;
D_Rs=0.568;
D_N=1.752;

k = 1.380649e-23;        #      Boltzmann constant
Q = 1.6021766208E-19;     #     Elementary charge of electron
Tref=300.15; #  with default temperature 27°C
D_Vt = k * Tref / Q;


function amps = Idiode(volts)
    global D_Vt D_N D_Is
    amps = D_Is * (exp(volts/(D_N*D_Vt))-1.0);
endfunction


# equations
# i1*R1 + v23 = V0
# V23 = i2*(R2+R3)
# i1 = i2 + i3
# i3*R4 + i3*D_Rs + Vd = V23
# i3 = Idiode(Vd)

# unknowns
#   1 i1
#   2 i2
#   3 i3
#   4 V23
#   5 Vd



function y = f (x)
    global R1 R2 R3 R4 V0 D_Rs
    y = zeros (5, 1);
    y(1) = x(1)*R1+x(4)-V0;
    y(2) = x(4)-x(2)*(R2+R3);
    y(3) = x(2)+x(3)-x(1);
    y(4) = x(3)*(R4+D_Rs) + x(5) -x(4);
    y(5) = x(3) - Idiode(x(5));
endfunction



[x, fval, info] = fsolve (@f, [0.001, 0.001, 0.001, 4.0, 0.4]);


x