function J=calc_finite_difference_jacobian(x,func,S,tol);

% Function to compute the Jacobian of a nonlinear function using 
% finite differences. 
% usage: J=calc_finite_difference_jacobian(x,func,[cols]);
%   x: (vector) value of independent variable
%   func: function handle for the nonlinear function. func(x) should 
%         return the nonlinear function as a column vector.
%   optional vector of column indices to limit computation to only certain columns
%   ouput J = d(func(x))/dx
% Notes: The Jacobian is computed one column at a time by calculating the 
% directional derivative J*e_j, where e_j is the unit vector in the j-th 
% direction. 
% Use this function if the Jacobian is dense or if it is sparse but you know nothing 
% about its sparsity pattern. If you do know the sparsity, use graph coloring 
% to compute the Jacobian more efficiently.
% This function uses the routine 'dirder.m' written by C. T. Kelley at NCSU. 
% The most recent version of this code appears to be as a nested function in:
% http://www4.ncsu.edu/~ctk/newton/SOLVERS/nsoli.m. You can extract it from 
% this file and save it as a separate m-file. Or you can email me for the 
% code.

% Samar Khatiwala (spk@ldeo.columbia.edu)

f0=func(x);

n=length(x);

useSparsity=0; % default
if nargin<3 % compute all columns
  cols=[1:n];
else
  if isempty(S) % assume want all columns of Jacobian
    cols=[1:n];
  elseif ~issparse(S)
    cols=S; % want some columns of Jacobian
  else
    useSparsity=1;
    if nargin<4
      tol=0;
    end
  end
end

if ~useSparsity
  m=length(cols);
  J=zeros(n,m); % return a dense matrix
  for j=1:m % loop over each column
    e=zeros(n,1); %  unit vector
    e(cols(j))=1;
    J(:,j)=dirder(x,e,func,f0); % (dfunc/dx)*e
  end
else
  g=colgroup(S);
  nr=length(unique(g));
  J=S;
  for ir=1:nr
    gr=find(g==ir);
    e=zeros(n,1);
    e(gr)=1;
    y=dirder(x,e,func,f0)+tol; % (dfunc/dx)*e
    for k=1:length(gr)
      irows=find(S(:,gr(k)));
      J(irows,gr(k))=y(irows);
    end
  end
end