Values for $a$ and $b$ as polynomials in $x,y$, when they exist, correspond to factorization of $\frac{x^p+y^2}{x+y}$ over the imaginary field $K_p:=\mathbb Q[\sqrt{-p}]$ of the form
$$\frac{x^p+y^2}{x+y} = (a(x,y)+\sqrt{-p}b(x,y))\cdot (a(x,y)-\sqrt{-p}b(x,y)).$$
[This Sage code](https://sagecell.sagemath.org/?z=eJyVkctqwzAQRfeB_MN0USK1bmynO1NnmX27TZ0w0QsVWxKKDIaQf-8Y4z4otHQ2EjNX9x5GUmkwKh3xxAKvlgugelk_DVuo4blHGTFZsbOqlewh8H0zKXY01SiSjwzYcAj3Jc9zNtABfFJYDc4nwLZlYNYdhqPwSmsrrHLpzIR3b73BpPhNbcBHkLebui5A09VkEqyjjJlnrKhSH6mJ7VlN3UgMIXpJ_gcm83zDv78mgl-DawrmMHvdfZrJ__lsv9hMkPGnvMXuJBFEBWJfNDz7S1M2tMfJdLkYccIIE6Lt1DGiM4o9ZmVRzBuigUtMry7hWr2my8eHXlf8HYF5kA8=&lang=sage&interacts=eJyLjgUAARUAuQ==) factors $\frac{x^p+1}{x+1}$ over $K_p$ and tests if the result has the required form, in which case the (univariate) polynomials $a(x,1)$ and $b(x,1)$ are reported. As an example, it reports $a$ and $b$ for all primes below $100$. Each prime $\equiv 3\pmod4$ in that range happens to have a solution.