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 this fact?  
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