p = 3
m = 1
F = GF(p)
Rx.<x> = PolynomialRing(F)
f = x^2 + 1
C_super = superelliptic(f, m)

Rxy.<x, y> = PolynomialRing(F, 2)
f1 = superelliptic_function(C_super, x^2)
f2 = superelliptic_function(C_super, x^4)
AS = as_cover(C_super, [f1, f2], prec=1000)
RxyzQ, Rxyz, x, y, z = AS.fct_field
f_z = as_function(AS, z[0]*z[1]*y)
df_z = f_z.diffn()
print(df_z.form == -z[0]*z[1]*x + y*z[1]*x - y*z[0]*x^3)