H0=-y*D[f,y]+x0*D[f,x0];
H1=-x0*D[f,x0]+x1*D[f,x1];
H2=-2*f-y*D[f,y]-x0*D[f,x0]-2*x1*D[f,x1];
a1=-x0*D[f,x1];
b0=y*D[f,x0];
b1=y*D[f,x1];
g0=x0*I*f;
g1=x1*I*f;
om=y*I*f;
am1=x1*D[f,x0];
bm0=-x0*D[f,y];
bm1=-x1*D[f,y];
gm0=-y*D[D[I*f,x0],y]-x0*D[D[I*f,x0],x0]-2*D[I*f,x0]-x1*D[D[I*f,x1],x0];
gm1=-y*D[D[I*f,x1],y]-x0*D[D[I*f,x1],x0]-x1*D[D[I*f,x1],x1]-2*D[I*f,x1];
mom=-y*D[D[I*f,y],y]-2*D[I*f,y]-x0*D[D[I*f,x0],y]-x1*D[D[I*f,x1],y];