%This code approximates the first derivatives of multivariate valued
%function. It requires as inputs the function handle and the points for
%which the derivative is to be approximated at. Written by Philip Shaw,
%Fordham University, 2022.
function[DCD]=CDJac(fx,xbar)
jk=length(xbar); %find the dimension of x
hstar=eps^(1/3); %choose a value of h based on Heer and Maussner machine epsilon
e=zeros(1,jk); %1 x jk vector of zeros
%loop for the numerical derivatives
for j=1:jk
e(j)=1; %replace the jth entry to 1 in the zero vector
fxbarph=fx([xbar+e.*hstar]); %function evaluated at point xbar plus h
fxbarmh=fx([xbar-e.*hstar]); %function evaluated at point xbar minus h
DCD(:,j)=(fxbarph-fxbarmh)./(2*hstar); %central difference derivative
e=zeros(1,jk); %create the ej row vector of zeros
end