function [I dIdU dIdC] = fni(U,C,param) tol = 1e-10; maxit = 10; delta = param.delta; E = param.E; R = param.R; T = param.T; F = param.F; alpha = param.alpha; i0 = param.i0; Cref = param.Cref; Rm = param.Rm; imax = C/delta; fail = 0; itspecial = 0; residual = 2*tol; I = 0.9*imax; nit = 0; maxit = 20; cutoff1 = 1e-3; cutoff2 = 0.999; res = E - R*T/alpha/F*(log(I/i0)- log((C-delta*I)/Cref))-Rm*I-U; while (residual > tol) && (fail ==0) residual = abs(res); % disp(sprintf('residual %d',residual)); dI = - R*T/alpha/F*(1/I + delta/(C-delta*I))-Rm; Iold = I; I = I - res/dI; if (I < cutoff1) && (itspecial ==0) % disp(sprintf('Special small current')); itspecial = 1; I = i0*C/Cref*exp(alpha*F/R/T*(E-U)); end if (I > cutoff2*imax) && (itspecial ==0) % disp(sprintf('Special large current')); itspecial = 1; temp = exp(alpha*F/R/T*(E-U-Rm*imax)); I = imax*(1-Cref/C*imax/i0/temp); end if I<= 0 % disp(sprintf('negative current correction')); I = Iold/2; end if I >=imax % disp(sprintf('more than maximum current correction')); I = Iold + (imax-Iold)/2; end nit = nit+1; if nit > maxit fail = 1; end if fail == 0 res = E - R*T/alpha/F*(log(I/i0)- log((C-delta*I)/Cref))-Rm*I-U; end end if fail ==1 disp(sprintf('fni iterations failed')) disp(sprintf('paused - hit control c')) pause end dIdU = 1/dI; dIdC = 1/(C/I+Rm*alpha*F/R/T*(C-delta*I));