function J=jacobian(v0)
  eps=1e-6; y0=shfun(v0);
  J=[ (shfun(v0+[eps;0])-y0)/eps, (shfun(v0+[0;eps])-y0)/eps ];
end
