Take $\Omega$ to be a bounded domain in $N$ dimensional Euclidean space with smooth boundary and we assume $\Omega$ contains the origin. I am interested is the following equation

$$ \Delta \phi(x) + \gamma \sum_{i,j=1}^N \frac{x_i x_j}{|x|^2} \phi_{x_i x_j}(x) = f(x)$$ in $\Omega$ with $ \phi=0$ on $ \partial \Omega$. Here $\gamma$ is such that the operator is elliptic. I would like to get some estimate of the form $$ \sup_{\Omega}| \nabla \phi(x)| \le C \sup_\Omega |f(x)|.$$ I guess there might be some issues at the origin. I have tried cutting the domain off and trying to get estimates independent of the domain cut off; but I was not successful. I guess part of my problem is that I don't truly understand the estimates I use; but normally I can perturb off the Lapclian and not really need to understand the main ideas.

QUESTION. Is this gradient estimate true and are there any standard methods to try and prove.
The method I tried was the $C^{1,\alpha}$ blow up analysis method from Leon Simon; but I couldn't get it to work.

thanks