function [psiN,dpsiN_dx] = eval_psiN(x,~,~,model,gauss)
% function [psi] = eval_psi(x,[],[],model,whichgrid)
%#codegen

psi = eval_psi(x,[],[],model,gauss);
if ~gauss
    psia = psi(1); psib=psi(end);
else
    psirho = eval_psi(x,[],[],model,false);
    psia = psirho(1); psib=psirho(end);
end
psiN = (psi-psia)/(psib-psia);

if ~gauss
    psiN(1) = 0;
end

if nargout>1
    error('to do')
    dpsiN_dx = zeros(numel(psi),numel(x));
    if ~gauss
        dpsiN_dx(:,model.psi.xind) = model.psi.Lam;
    else
        dpsiN_dx(:,model.psi.xind) = model.psi.Lamgauss;
    end
end


return
