function [J] = jacfd(F, x, workers)
% function [J] = jacfd(F, x[, workers])
% Finite difference jacobian evaluation
% F: Function handle which takes x as argument
% x: F(x) point around which to evaluate Jacobian
% workers: number of workers used. Deafult:0, no parallel computation
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

if nargin<3
  workers = 0; % don't use parallel computation by default
end
% Compute the Jacobian matrix via central finite differences

nn = numel(x);
J = zeros(nn,nn);
doplot = false;
epsval = sqrt(eps);
II = eye(nn);

if workers>1
  parpool(1,workers)
  parfor jj=1:nn
    dx = epsval*II(:,jj);
    % central differences
    J(:,jj) = ( F( x + dx,doplot ) - F(x - dx,doplot ) )/2/epsval;
  end
else
  % no parallel  computation
  for jj=1:nn
    dx = epsval*II(:,jj);
    % central differences
    J(:,jj) = ( F( x + dx,doplot ) - F(x - dx,doplot ) )/2/epsval;
  end
end

assert(~any(isnan(J(:))), 'NaNs in Jacobian');
end