function [dce, dci] = GRADFUN(x)
% finds finite difference gradients corresponding to FUN.m
global nvars ne ni

% dce(i,j) = derivative of constraint j w.r.t variable xi

if ne > 0, 
	dce = zeros(nvars,ne); 
	ce1 = zeros(ne,1);
	ce2 = zeros(ne,1);
end
dci = zeros(nvars,ni);

h = 1e-4;
if ni> 0,
ci1 = zeros(ni,1);
ci2 = zeros(ni,1);
end

for i = 1:nvars
	xplus = x;
	xplus(i) = x(i) + h;
	xminus = x;
	xminus(i) = x(i) - h;

	[f, ce1, ci1] = FUN(xminus);
	[f, ce2, ci2] = FUN(xplus);

% 	using central differences

	if ce1 == [],
	   dce = [];
	else
	   dce(i,1:ne) = (ce2' - ce1')/h ;
	end
	if ci1 == [],
	   dci = [];
	else
	   dci(i,1:ni) = (ci2' - ci1')/h ; % ci never empty
	end
    
end
