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));