# A better optimized version of DF
function M=DFopt(x,y)
      t2 = x^2;
      t3 = t2*x;
      t4 = 3*t3;
      t5 = y^2;
      t6 = x*t5;
      t7 = 9*t6;
      t9 = t2*y;
      t10 = 3*t9;
      t12 = t5*y;
      t14 = 3+t2+t5;
      t15 = t14^2;
      t17 = 1/t15/t14;
      d2zdx2 = 2*(-27*x+t4-t7+3*t2-t10-3-t5+3*y+t12)*t17;
      t19 = 9*y;
      d2zdy2 = -2*(-t19-t10+t12-3*t5-t7+3+t2+9*x+t4)*t17;
      d2zdxdy = 2*(-t19+9*t9-3*t12+3*x+t3-3*t6+4*x*y)*t17;
      M(1,1)=d2zdx2;
      M(2,2)=d2zdy2;
      M(2,1)=d2zdxdy;
      M(1,2)=d2zdxdy;
end
