WITH
  (* Derivatives *)
  f = too,
  fx = (tpo - tmo)/2,
  fy = (top - tom)/2,
  fxy = (-tmp+tpp+tmm-tpm)/4,
  
  sd = tmp+tpp+tmm+tpm,
  so = tmo+tpo+tom+top+too,
  
  sk = (sd-so-so)/6,

  fyy = sk + (top + tom - too),
  fxx = sk + (tpo + tmo - too),
  
  (* Compute E, dEdf, dEdfx, ..., dEdfyy *)
  
  dEdsk  += dEdfxx;
  dEdtpo += dEdfxx;
  dEdtmo += dEdfxx;
  dEdtoo += - dEdfxx;
  
  dEdsk  += dEdfyy;
  dEdtop += dEdfyy;
  dEdtom += dEdfyy;
  dEdtoo += - dEdfyy;
  
  dEdsd += dEdsk/6;
  dEdso += - dEdsk/3;
  
  dEdtmo += dEdso;
  dEdtpo += dEdso;
  dEdtom += dEdso;
  dEdtop += dEdso;
  dEdtoo += dEdso;
  
  dEdtmm += dEdsd;
  dEdtmp += dEdsd;
  dEdtpm += dEdsd;
  dEdtpp += dEdsd;
  
  dEdtmp += - dEdfxy/4,
  dEdtpp += dEdfxy/4,
  dEdtmm += dEdfxy/4,
  dEdtpm += - dEdfxy/4;
  
  dEdtop += dEdfy/2,
  dEdtom += - dEdfy/2;
  dEdtpo += dEdfx/2;
  dEdtmo += - dEdfx/2;
  
  dEdtoo += dEdf
  
  
DO
  ...
END

WITH
  gx = fx*fx,
  gy = fy*fy,
  gm = 2*fx*fy
  
  p0 = f,
  p1 = gx + gy
  p2 = fxx + fyy,
  p3 = fxx*gy - fxy*gm + fyy*gx,
  p4 = fxy*fxy - fxx*fyy,
  
  (* Now compute E, dEdp *)
  
  dEdfyy += - dEdp[4]*fxx;
  dEdfxx += - dEdp[4]*fyy;
  dEdfxy += dEdp[4]*fxy*2;
  
  dEdfxx += dEdp[3]*gy;
  dEdgy  += dEdp[3]*fxx;
  dEdfxy += - dEdp[3]*gm;
  dEdgm  += - dEdp[3]*fxy;
  dEdfyy += dEdp[3]*gx;
  dEdgx  += dEdp[3]*fyy;
  
  dEdfxx += dEdp[2];
  dEdfyy += dEdp[2];
  
  dEdgx += dEdp[1];
  dEdgy += dEdp[1];

  dEdf  += dEdp[0];
  
  dEdfx += dEdgm * 2 * fy;
  dEdfy += dEdgm * 2 * fx
  dEdfy += dEdgy * 2 * fy;
  dEdfx += dEdgx * 2 * fx;
  
END;

(* Program to Maple *)
  
(* InvariantsFromDerivatives: *)

f := v[0];
fx := v[1];
fy := v[2];
fxx := v[3];
fxy := v[4];
fyy := v[5];

gx = fx*fx;
gy = fy*fy;
gm = 2.0*fx*fy;

p[0] := f;
p[1] := gx + gy;
p[2] := fxx + fyy;
p[3] := fxx*gy - fxy*gm + fyy*gx;
p[4] := fxy*fxy - fxx*fyy;

(* InvariantPartialsFromDerivativePartials: *)

gx := v[1]*v[1];
gy := v[2]*v[2];
gm := 2*v[1]*v[2];
dEdgm := - dEdp[3]*v[4];
dEdgx := dEdp[1] + dEdp[3]*v[5];
dEdgy := dEdp[1] + dEdp[3]*v[3];
dEdv[0] := dEdp[0];
dEdv[1] := 2*(dEdgm * v[2] + dEdgx * v[1]);
dEdv[2] := 2*(dEdgm * v[1] + dEdgy * v[2]);
dEdv[3] := dEdp[2] + dEdp[3]*gy - dEdp[4]*v[5];
dEdv[4] := - dEdp[3]*gm + 2*dEdp[4]*v[4];
dEdv[5] := dEdp[2] + dEdp[3]*gx - dEdp[4]*v[3];

   (* Program to Maple *)
   
dEdf := dEdv[0];
dEdfx := dEdv[1];
dEdfy := dEdv[2];
dEdfxx := dEdv[3];
dEdfxy := dEdv[4];
dEdfyy := dEdv[5];
dEdsk := dEdfxx + dEdfyy;
dEdsd := dEdsk/6;
dEdso := - dEdsd - dEdsd;
HdEdfx := dEdfx/2;
HdEdfy := dEdfy/2;
QdEdfxy := dEdfxy/4;
dEdt[-1,-1] := dEdsd + QdEdfxy;
dEdt[+1,-1] := dEdsd - QdEdfxy;
dEdt[-1,+1] := dEdsd - QdEdfxy;
dEdt[+1,+1] := dEdsd + QdEdfxy;
dEdt[00,-1] := dEdfxx - HdEdfx + dEdso;
dEdt[00,+1] := dEdfxx + HdEdfx + dEdso;
dEdt[-1,00] := dEdfyy - HdEdfy + dEdso;
dEdt[+1,00] := dEdfyy + HdEdfy + dEdso;
dEdt[00,00] := dEdf - dEdfxx - dEdfyy + dEdso;