%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 
    
    
