function J = numjac(mvec)
% usage: J = numjac(mvec)
% description: computes the numerical jacobian of a vector
% at the argument mvec.  Here the function name myfcn is 
% hardwired into this routine, as is the stepsize h. So
% user must adjust these before using this routine.

% we want to avoid using function name as extra argument
% so define myfcn in a separate function file and change names
% here if you use a different name for the function.
h = 1e-7; % put stepsize of your choice here
fcnmvec = myfcn(mvec);
m = length(fcnmvec);
n = length(mvec);
J=zeros(m,n);
Id = eye(n);
for ii=1:n
  J(:,ii)= (myfcn(mvec + h*Id(:,ii)) - fcnmvec)/h;
end
