current density in metallic homogeneous sphere You have a metallic homogeneous sphere. At the left side of the sphere electric current goes inside and at the right the same amount of current goes out of the sphere. The question is: What is the current density j(x) inside the sphere?

My reasoning:
Let's say that sphere has a radius R and its center is at $x = 0$. Origin of the system is at $(x = 0, y = 0)$. Current density is: $$j = \frac{I}{S}.$$ In the case of sphere the surface through which current if flowing is a function of $x$ which means that the current density is also function of $x$. I think the current density for this problem is then: $$j(x) = \frac{I}{\pi*(R^2 - x^2)}.$$
This solution however gives singularities at $j(R)$ and $j(-R)$ so I am a bit sceptical.
 A: The question boils down to the solution of the Laplace equation for the electrostatic potential $\phi$ in the sphere with the given boundary conditions.
First, you have to recognize that from the current continuity equation $$div{ \vec j}=-\frac{\partial \rho}{\partial t}=0$$ together with the current equation $$\vec j=\sigma \vec E$$ and $\vec E=-grad\phi$  you obtain the Laplace equation for the potential $\phi$ $$∆ {\phi}=0$$ with the boundary conditions normal current density zero, i.e., $$\vec E_n =-\vec n ·grad{\phi}=0$$ on the surface of the sphere with the exception of the current entrance and exit points and the potentials at the current entrance and exit points $\phi(in)=V$ and $\phi(out)=0$, respectively,  where  $V=I·R$ is the applied voltage following from the given current $I$ and the calculated resistance $R$ of the sphere. 
A: Your singularities are perfectly reasonable in my opinion and they come from a very common idealisation in electrodynamics: the assumption of "thin" wires with no diameter, of "point" charges with no radius and surfaces without thickness. In any of these cases, you will encounter a similar singularity to the one you have here. Try e.g. to think about a sphere with charge $q$ and radius $r$, therefore volume $\frac{4}{3}\pi r^3$. Its charge density is just $\varrho=\frac{q}{V}$. In the limit $r \to 0$ (a point charge) $V\to0$ and the charge density diverges because we want to put charge into a space without volume so to speak. There are two ways to deal with this. The usual one is to use Dirac/delta distributions which mathematically provide us with an object that does what we want: its volume integral equals the charge and it's "zero everywhere but that point" (although it's not really a function ...). The more complicated way is to give everything a finite volume (e.g. the wires on both sides of your sphere that provide the current into the sphere). Then you avoid infinities, but things often become very complex.
For the actual form of the current density, I think Elliot Yu's answer is very thoroughly thought through.
A: First of all, to directly answer your questions about your own approach, it is sensible that you would get singularities at $R$ and $−R$. Those are the points where $j$ enters and leaves the sphere, and the current density is concentrated at a point, so you would get a delta function-like distribution. However your assumption that the current density is only along the $x$-axis and that it is uniform on the cross-sectional surface both seem unreasonable.
Here I will present a brief outline of how I would go about solving this problem. I will skip over some computational details and the final solution for all of the components of the current density. Any edits to fill in details are welcome. I will make use  of Green's functions and the representation formula for a sphere.[1]
First due to Ohm's law, we have $\mathbf j = \sigma \mathbf E$. Under DC conditions, we have $\nabla\times \mathbf E = 0$. We can apply the curl operator on this equation again, and get $\nabla(\nabla\cdot \mathbf E) - \nabla^2 \mathbf E = 0$. By Gauss's law the first term drops out. Finally we end up with the formulation of the problem that $\nabla^2 \mathbf j = 0$. Note that this is a vector Laplace's equation, and due to the symmetry of the problem, we would like to choose the cylindrical coordinates with the $z$ axis parallel to the direction of the entering and exiting currents, so the equation actually takes a more complicated form. Using symmetry arguments, we can also see that $\mathbf j$ has no $\phi$ components, and the $z$ and $\rho$ components have no $\phi$ dependence. Thus the equation simplifies to the following two,
\begin{gather}
\nabla^2 j_z = 0 \>,\\
\nabla^2 j_\rho - \frac{j_\rho}{\rho^2} = 0 \>.
\end{gather}
We need to also find the boundary condition. Since we are considering a DC current, the electric fields are confined to the conductor, and the outside fields are zero, which then means that on the boundary $\mathbf j\sim \mathbf E = 0$, except for the entry and exit points. Since at those points we have a current $I$, the boundary condition should be $\mathbf j|_{\text{boundary}}(\mathbf r) = I(\delta^2(z - R) + \delta^2( z +R))\hat{z}$. Here $z = R$ is the position of the exit point. Note that here the delta functions are only two dimensional, because of the definition of $\mathbf j$ as the current passing through unit area, $I = \int_S \mathbf j \cdot d\mathbf S$.
With all of these under our belt, we can then use the representation formula for Laplace's equation with Dirichlet boundary condition to solve the $j_z$ equation, (here written in spherical coordinates for simplicity of notations)
$$
j_z(\mathbf r) = \frac{R^2-r^2}{4\pi R}\int_{\text{boundary}} \frac{g(\mathbf r')}{|\mathbf r - \mathbf r'|^3}dS' \>.
$$
If we integrate this out we should get
$$
j_z (\rho, z) =\frac{R^2-(\rho^2 + z^2)}{4\pi R}\left(\frac{I}{(\rho^2+(R-z)^2)^{3/2}}+\frac{I}{(\rho^2+(R+z)^2)^{3/2}}\right) \>.
$$
Now that we know $j_z$, there is actually no need to also solve the $j_\rho$ equation, which is arguably more difficult. Due to charge conservation, and the fact that there is no net charge inside the conductor, we can use continuity equation to find $j_\rho$, that is,
$$
\frac{1}{\rho}\frac{\partial(\rho j_\rho)}{\partial \rho} + \frac{\partial j_z}{\partial z} = 0 \>.
$$
Note: Somehow the result does not seem quite right. The integral $\int_S \mathbf j \cdot d\mathbf S$ does not equal $I$ at $z=0$. The issue might have come from the normalization of the delta functions and/or the way the integral was performed on the delta functions. The general principle should still hold, but I shall try to correct my calculations.
[1] L. C. Evans, Partial Differential Equations (American Mathematical Society, Providence, RI, 2010), pp. 33-41.
