%This code approximates the first derivative of a 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, 2018.
function[DCD]=CDJac(fx,xbar)
jk=length(xbar); %find the dimension of x
hstar=eps^(1/3); %choose value of h based upon Heer and Maussner machine epsilon
e=zeros(1,jk); %1 x j row vector of zeros
for j=1:length(xbar)
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 estimate for the derivative
e=zeros(1,jk); %create the ej row vector of zeros
end