function [res, LY] = liuF(x, L, LX, x0, doplot) 
% Operator for standard LIUQE with JFNK
% Warning dz stabilization is kept only for 1:1 comparison with LIUQE standard but can be removed
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

indIe  =               (1:L.G.na+L.G.nu);
indIyD = indIe (end) + (1:L.ny*L.nD);
inddz  = indIyD(end) + (1:L.ndz);

% Extract the inputs
Ie  = x(indIe );
IyD = x(indIyD);

if L.P.stabz
  dz = x(inddz);
else
  dz = 0;
end

IyD = reshape(IyD, L.nzy, L.nry , L.nD);

LY = liut(L,LX,'rst',false,'IyD',IyD, 'Ie', Ie, 'dz', dz);

if isempty(LY)
  res = NaN;
else
  % Compute residual
  resIe = (Ie - [LY.Ia;LY.Iu])/max(abs(x0(indIe)));
  resIyD = reshape(IyD - LY.IyD, L.ny*L.nD, 1)/max(abs(x0(indIyD)));
  
  res = [ resIe; resIyD];
  if L.P.stabz
    resdz = (dz - LY.dz);
    res = [res; resdz];
  end
end
  

if doplot && all(~isnan(res))
  clf
  % Debugging plot if requested
  meqplotliu(L,LX,LY, gcf)
  drawnow
end

end

