function y=F(x)
PincdBm=25;
Pinc=10^(PincdBm/10-3); // conversion to W
Io=0;
Is=54e-9;
Rs=7;
L=38.5;
n=1.12;
Rg=50;
RL=500;
// besseli is the modified bessel function
y=besseli(0,L/n*sqrt(8*Rg*Pinc))-(1+Io/Is+x/(RL*Is))*exp((1+(Rg+Rs)/RL)*L/n*x+L/n*Rs*Io);
endfunction
err=1e-3;
x0=0.6735; // first guess
maxiter=5000; // max iterations
iter=0;
delta=1e-2;
while (abs(F(x0))>err)
iter=iter+1;
if (F(x0-delta)*F(x0+delta) < 0) then
delta=delta/1.5;
elseif ((F(x0+delta)-F(x0-delta))>0) then
x0=x0-delta;
delta=delta*2;
else
x0=x0+delta;
delta=delta*2;
end
if ((iter>=maxiter)&(F(x0-delta)*F(x0+delta) < 0)) then
break;
end
end
x0 // print the estimated root
delta // print the error range (the true zero is between x0-delta and x0+delta)