Given variables $x \in \mathbb{R}^{d_1}$, $y \in \mathbb{R}^{d_2}$ with $d_1 + d_2 = d$ and smooth functions $f_1, f_2 : \mathbb{R}^d \to \mathbb{R}$, I want to solve the equation $$g_1(x,y) = -(\nabla_x f_1)(x, y-\nabla_y f_2(x+g_1(x,y), y))$$ in $g_1 : \mathbb{R}^{d} \to \mathbb{R}^{d_1}$, where the RHS notation means the gradient of $f_1$ evaluated at the point $(x, y-\nabla_y f_2(x+g_1(x,y), y))$.
This seems intractable in general, but all I need is to approximate $g_1$ numerically at any given point $(x,y)$ for any given functions $f_1, f_2$. I am not familiar with differential equations/functions defined implicitly; I heard about finite difference methods but can anyone help me through how to use them here? In particular I'm looking to implement this in Python.