function [Y,inewton]=BE(Y,dt,Iinj,tol,inewtonmax,nskip)

  U=Y; 
  inewton=1;
  diffy=1;   
    
  iskip=0;
  while diffy > tol
     
    if mod(iskip,nskip) == 0
      jacobian=dG_BE(dt,U,Y,Iinj);
    end
    iskip=iskip+1;
    gbe=G_BE(dt,U,Y,Iinj)';
    dY=jacobian\gbe;
    U=U-dY';
    diffy=norm(dY);
    inewton=inewton+1;
%    disp([inewton(k),diffy]);
    if inewton > inewtonmax 
     disp(' Newton does not converge');
%     [t(k),u,dy,inewton(k)]
      break
    end
    
  end 
  
  Y=U;

end