restart;
z:=(1+3*x-y)/(3+x^2+y^2);

with(plots):
# Warning, the name changecoords has been redefined
contourplot(z,x=-3..3,y=-3..3);

dzdx:=diff(z,x);
dzdy:=diff(z,y);


with(codegen):
# Warning, the protected name MathML has been redefined and unprotected
fortran(dzdx);
#       t0 = 3/(3+x**2+y**2)-2*(1+3*x-y)/(3+x**2+y**2)**2*x
dzdxo:=[optimize(simplify(dzdx))];

fortran(dzdxo);
#       t1 = x**2
#       t3 = y**2
#       t10 = (3+t1+t3)**2
#       t13 = -(-9+3*t1-3*t3+2*x-2*x*y)/t10
fortran([optimize(simplify(dzdy))]);
#       t1 = x**2
#       t2 = y**2
#       t8 = (3+t1+t2)**2
#       t11 = -(3+t1-t2+2*y+6*x*y)/t8
d2zdx2:=diff(dzdx,x);
d2zdy2:=diff(dzdy,y);
d2zdxdy:=diff(dzdx,y);



dzdx;

evalf(subs({x=1,y=1},dzdx));

fortran([optimize(simplify(d2zdx2))]);
#       t2 = x**2
#       t5 = y**2
#       t14 = 3+t2+t5
#       t15 = t14**2
#       t19 = 2*(-27*x+3*t2*x-9*x*t5+3*t2-3*t2*y-3-t5+3*y+t5*y)/t15/t14
fortran([optimize(simplify(d2zdy2))]);
#       t2 = x**2
#       t5 = y**2
#       t14 = 3+t2+t5
#       t15 = t14**2
#       t20 = -2*(-9*y-3*t2*y+t5*y-3*t5-9*x*t5+3+t2+9*x+3*t2*x)/t15/t14
fortran([optimize(simplify(d2zdxdy))]);
#       t2 = x**2
#       t5 = y**2
#       t15 = 3+t2+t5
#       t16 = t15**2
#       t20 = 2*(-9*y+9*t2*y-3*t5*y+3*x+t2*x-3*x*t5+4*x*y)/t16/t15
fortran([optimize(['d2zdx2'=simplify(d2zdx2),'d2zdy2'=simplify(d2zdy2),'d2zdxdy'=simplify(d2zdxdy)])]);
#       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