I have a function $u(x,y)$ describing a certain temperature distribution that looks like
$u(x,y) = -\dfrac{q}{4\pi \lambda} [\ln((x-a)^2+(y-a)^2)-\ln((x+a)^2+(y-a)^2)+\ln((x+a)^2+(y+a)^2)-\ln((x-a)^2+(y+a)^2)], \quad \text{for} \:x>0,y>0$
and it is given that $u(0,y) = 0$, $u(x,0) =0$.
Now I am interested in finding the maximum points for the heat flux density $\textbf{j} =-\lambda\nabla u$ going out from the boundaries. So I thought that shouldn't be a problem, we just solve the equations $|\nabla_y u(x,y)|_{y=0} = 0$ and $|\nabla_x u(x,y)|_{x=0} = 0$ (this would give us the magnitudes of the normal derivatives at the boundaries, if you disagree then I am interested to hear why). But according to my calculations, both of these expressions lead to $0=0$ (which is probably obvious in some way and if not then it is obvious my calculations are inaccurate). So I ask you, the reader, for a push in the right direction.
you will get $$\frac{\partial f(x,y)}{\partial x}=-1/4\,{\frac {q}{\pi \,\lambda} \left( {\frac {2\,x-2\,a}{ \left( x-a \right) ^{2}+ \left( y-a \right) ^{2}}}-{\frac {2\,x+2\,a}{ \left( x+ a \right) ^{2}+ \left( y-a \right) ^{2}}}+{\frac {2\,x+2\,a}{ \left( x +a \right) ^{2}+ \left( y+a \right) ^{2}}}-{\frac {2\,x-2\,a}{ \left( x-a \right) ^{2}+ \left( y+a \right) ^{2}}} \right) } $$ and $$\frac{\partial f(x,y)}{\partial y}=-1/4\,{\frac {q}{\pi \,\lambda} \left( {\frac {2\,y-2\,a}{ \left( x-a \right) ^{2}+ \left( y-a \right) ^{2}}}-{\frac {2\,y-2\,a}{ \left( x+ a \right) ^{2}+ \left( y-a \right) ^{2}}}+{\frac {2\,y+2\,a}{ \left( x +a \right) ^{2}+ \left( y+a \right) ^{2}}}-{\frac {2\,y+2\,a}{ \left( x-a \right) ^{2}+ \left( y+a \right) ^{2}}} \right) }$$