H0=x0*D[f,y];
H1=-x0*D[f,x0]+x1*D[f,x1];
H2=-x1*D[f,x1]+x2*D[f,x2];
H3=-x2*D[f,x2]+x3*D[f,x3];
a1=-x0*D[f,x1];
a2=x1*D[f,x2];
a3=x2*D[f,x3];
a4=-x0*D[f,x2];
a5=x1*D[f,x3];
a6=-x0*D[f,x3];
b0=y*D[f,y];
b1=y*D[f,x1];
b2=y*D[f,x1];
b3=y*D[f,x1];
g0=x0*I*f;
g1=x1*I*f;
g2=x2*I*f;
g3=x3*I*f;
om=y*I*f;
gm0=-y*D[D[I*f,x0],y]-x0*D[D[I*f,x0],x0]-3*D[I*f,x0]-x1*D[D[I*f,x1],x0]-x2*D[D[I*f,x2],x0]-x3*D[D[I*f,x3],x0];
am1=x1*D[f,x0];
am2=-x2*D[f,x1];
am3=-x3*D[f,x2];
am4=x2*D[f,x0];
am5=-x3*D[f,x1];
am6=x3*D[f,x0];
bm0=-x0*D[f,y];
bm1=x1*D[f,y];
bm2=x2*D[f,y];
bm3=x3*D[f,y];
gm1=y*D[D[I*f,x1],y]+x0*D[D[I*f,x1],x0]+x1*D[D[I*f,x1],x1]+3*D[I*f,x1]+x2*D[D[I*f,x2],x1]+x3*D[D[I*f,x3],x1];
gm2=y*D[D[I*f,x2],y]+x0*D[D[I*f,x2],x0]+x1*D[D[I*f,x2],x1]+x2*D[D[I*f,x2],x2]+3*D[I*f,x2]+x3*D[D[I*f,x3],x2];
gm3=y*D[D[I*f,x3],y]+x0*D[D[I*f,x3],x0]+x1*D[D[I*f,x3],x1]+x2*D[D[I*f,x3],x2]+x3*D[D[I*f,x3],x3]+3*D[I*f,x3];
mom=-y*D[D[I*f,y],y]-3*D[I*f,y]-x0*D[D[I*f,x0],y]-x1*D[D[I*f,x1],y]-x2*D[D[I*f,x2],y]-x3*D[D[I*f,x3],y];