function [LY] = liutnwt(L,LX)

%% Set default parameters for the JFNK solver
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
Pjfnk.mkryl = 200;
Pjfnk.tolF =  1e-10;
Pjfnk.robust_on =  1;
Pjfnk.kmax = 100;
Pjfnk.debug = L.P.debug;

%% Initialize restarting flag
solverinfo.isconverged = 0;
ii = 1; % Counter for converged equilibria
LYt = []; % Temporary LY as output of the solver to be checked


%% Define L structures
if ~isfield(L.P,'argout'), L.P.argout = {};end
% L structure for liut during for re-initialization
P_restart  = L.P;
P_restart.itert = 3; % Do not take the iterated quantities from input but restart from FE
P_restart.iters = 0; % Use the backslash linear solver
P_restart.argout{end+1} = 'IyD';
L_restart = liuc(P_restart,L.G);

% L structure for liut during Newton steps
P_iteration = L.P;
P_iteration.debug = 0; % Do not use the liut debugging tools ( otherwise computed at every krylov direction and every Newton step)
P_iteration.itert = 1; % Initialize from iterated quantities provided externally
P_iteration.iters = 0; % Use the backslash linear solver
P_iteration.argout{end+1} = 'IyD';
L_iteration = liuc(P_iteration,L.G);


%% Time loop
for tt=1:numel(LX.t)
 % Initial condition from FE
 if tt == 1 || solverinfo.isconverged == 0
  if abs(LX.Ip(tt))< L.P.Ipmin
   continue
  end
  % Restart from FE and itera Picard iterations
  LYt = liut(L_restart,meqxk(LX,tt));
 end
 
 if ~isempty(LYt) % Start running only when initialization is available
  % Update the initial condition for the solver
  x0 = [[LYt.Ia;LYt.Iu]; LYt.IyD(:); LYt.dz(1:L.ndz)];
  
  % Define NL to be solved
  F = @(x,doplot) liuF(x,L_iteration,meqxk(LX,tt),x0, doplot);
  
  % Solve the problem with JFNK
  [x1, solverinfo] = solveF(F,x0, 'debug', Pjfnk.debug,...
   'mkryl', Pjfnk.mkryl, 'tolF', Pjfnk.tolF, 'robust_on',...
   Pjfnk.robust_on,'kmax', Pjfnk.kmax );
  
  %  Compute post-processing
  [res,LYt] = F(x1,0);
  
  % Solution sanity check
  if solverinfo.isconverged
   % Store results
   LY(ii) = LYt;
   ii=ii +1;
  end
 end
 
end

if ii>1 % At least 1 solution was found
 % Pack results in usual LIUQE form
 LY = meqlpack(LY,{'ppQg','TTpQg','ag'});
 LY.shot = LX.shot;
else
 LY = [];
end

end

 
