BL_f.m

function f=BL_f(u,M)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%           water fractional flow function f(u)   --   BL             %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% function f=BL_f(u)


f = zeros(size(u));

for k=1:length(u)
    if (u(k)>=0) && (u(k)<=1)
        f(k) = u(k)^2/(u(k)^2+M*(1-u(k))^2);
    elseif u(k)<0
        f(k) = 0;
    else
        f(k) = 1;
    end
end

end