function [ l, u ] = diffjac ( x, f, f0 )

%*****************************************************************************80
%
%% diffjac() computes and factors the jacobian of a function.
%
%  Discussion:
%
%    A forward difference estimate is used.
%
%    The routine dirder() is needed to compute the directional derivative
%    that defines each column of the jacobian estimate.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    13 April 2021
%
%  Author:
%
%    Tim Kelley
%
%  Reference:
%
%    Tim Kelley,
%    Iterative Methods for Linear and Nonlinear Equations,
%    SIAM, 2004,
%    ISBN: 0898713528,
%    LC: QA297.8.K45.
%
%  Input:
%
%    real X(N): the point at which the jacobian is estimated.
%
%    function fx = F ( x ): MATLAB code that evaluates the function.
%
%    real F0(N): the value of F(X).
%
%  Output:
%
%    real L(N,N), U(N,N): the LU factorization of the (approximate) jacobian.
%
  x = x(:);

  n = size ( x, 1 );

  for j = 1 : n
    zz = zeros ( n, 1 );
    zz(j) = 1;
    jac(:,j) = dirder ( x, zz, f, f0 );
  end

  [ l, u ] = lu ( jac );

  return
end

