function H = mpecd2obj(x,y)
% H = mpecd2obj(x,y)
% dense finite difference approx to Hessian of theta

[theta,gx,gy,domerr] = mpecfunc(4,x,y);
lintheta = mpecfunc(7,x,y);

z = [x; y];
n = length(z);
nx = length(x);
delta = sqrt((eps/2)^(0.8));
H = sparse(n,n);
for j=1:n
  if (lintheta(j) == 0)
    zj = z(j);
    dz = delta*(1.0+abs(zj));
    z(j) = zj + dz;
    [theta,gxU,gyU,domerr] = mpecfunc(4,z(1:nx),z(nx+1:n));
    H(:,j) = sparse([gxU-gx; gyU-gy] / dz);
    z(j) = zj;
  end
end;

return;

% the following code does central differencing as opposed to 
% finite differencing
z = [x; y];
n = length(z);
nx = length(x);
delta = (eps/2)^(1/3);
for j=1:n
  zj = z(j);
  dz = delta*(1.0+abs(zj));
  z(j) = zj + dz;
  [theta,gxU,gyU] = mpecfunc(4,z(1:nx),z(nx+1:n));
  z(j) = zj - dz;
  dz = dz + dz;
  [theta,gx,gy] = mpecfunc(4,z(1:nx),z(nx+1:n));
  h = [gxU-gx; gyU-gy] / dz;
  z(j) = zj;
  H = [H h];
end;

