# Recursive Function in Multiple Variables

I can define a recursive function in one variable, e.g. defining the factorial function $$f(n)= n!$$ using Mathematica. But how do you do this for a function which depends not only on $$n$$ but on $$x$$ and $$y$$. Ultimately, I would like to define the Division Polynomials, which depend on the $$n$$ (both an odd and even case) and then inputs $$x,y$$, which the function will be evaluated at later.

(*" This code follows the Division Polynomial Wikipedia article "*)
ClearAll[s, d, psi, phi, omega, P, x, y, X, Y, Z, A, B];
s[0] = d[0] = 0; d[1] = d[2] = 1; d[3] = Y; d[4] = Z;
s[n_ /; n < 0] := -s[-n];
d[n_ /; n < 0] := -d[-n];
s[n_Integer] := If[OddQ[n], 1, X] d[n];
d[n_ /; OddQ[n]] := d[n] = With[{m = (n - 1)/2},
If[OddQ[m], 1, X^4] d[m + 2] d[m]^3 -
If[OddQ[m], X^4, 1] d[m - 1] d[m + 1]^3 // Factor];
(*" This version involves **no** division "*)
d[n_ /; EvenQ[n]] := d[n] = With[{m = n/2},
d[m] (d[m + 2] d[m - 1]^2 - d[m - 2] d[m + 1]^2) // Factor];
(*" This is the division polynomial and related functions "*)
psi[n_Integer, x_:x, y_:y] := psi[n, x, y] = s[n] /. {X -> 2 y,
Y -> 3 x^4 + 6 A x^2 + 12 B x - A^2,
Z -> 2 (x^6 + 5 A x^4 + 20 B x^3 - 5 A^2 x^2 - 4 A B x - 8 B^2 - A^3)};
phi[n_Integer, x_:x, y_:y] := phi[n, x, y] =
x psi[n, x, y]^2 - psi[n + 1, x, y] psi[n - 1, x, y];
omega[n_Integer, x_:x, y_:y] := omega[n, x, y] =
(psi[n + 2, x, y] psi[n - 1, x, y]^2 -
psi[n - 2, x, y] psi[n + 1, x, y]^2)/(4 y) // Factor;
P[n_Integer, x_:x, y_:y] := p[n, x, y] =
{phi[n, x, y]/psi[n, x, y]^2, omega[n, x, y]/psi[n, x, y]^3};


Test it with, for example, P[1] == {x, y} or P[1, x, y] == {x, y}.