I am trying to solve for the function $p(x)$ which obeys the following: \begin{align*} p\left(x\right) = \left(e^{x} - \frac{1}{2} z \left(p\left(x\right)\right)^2\right)^\gamma \phi(x), \tag{1} \end{align*} where $\phi(x)$ satisfies \begin{align*} e^x \left[e^{x} - \frac{1}{2} z \left(p\left(x\right)\right)^2\right]^{-\gamma} = \left[\lambda + (1-\gamma)\delta + \gamma z p(x)\right]\phi(x) - \mathcal A\left(\phi(x) \right), \tag{2} \end{align*} and $\mathcal A\left(\cdot \right)$ is the infinitesimal generator specified as: \begin{align*} \mathcal A\left(g(x) \right) = - \theta_x \left(x - \bar x\right)g'(x) + \frac{1}{2} \sigma_x^2 g''(x) \tag{3} \end{align*}
I use Python to approximate $p(x)$ with a polynomial and use the optimizer to find the appropriate coefficients, but the optimizer is highly sensitive even on the small range of $x \in [-4.8, -4.2]$ I am looking at, and it converges (or diverges..) to extremely different values depending on the initial guess. I was hoping someone could suggest a proper initial guess for the coefficients as well as an appropriate degree of the polynomial.
An alternative solution is also welcome. I tried to rewrite the problem as a system of equations \begin{align*} p(x) = f(p(x), \phi(x)) = \left(e^{x} - \frac{1}{2} z \left(p\left(x\right)\right)^2\right)^\gamma \phi(x) \tag{4} \end{align*} and \begin{align*} \phi(x) = h(p(x), \phi(x)) = e^{-x}p(x)\bigg(\left[\lambda + (1-\gamma)\delta + \gamma z p(x)\right]\phi(x) - \mathcal A\left(\phi(x) \right)\bigg), \tag{5} \end{align*} and then iterate on $f(\cdot)$ and $h(\cdot)$, but neither of them are contractions, so this method doesn't seem to work.
Thanks!