This is basically a problem of recursive counting. Start with the uncoupled basis state, i.e. the set of states of the form
$\vert \ell m_\ell \rangle \vert s m_s\rangle$. There are clearly $(2\ell+1)(2s+1)$ of these, and the job is to reorganize them.
The key counting result is based on the observation that
$\vert \ell m_\ell\rangle\vert sm_s\rangle$ is an eigenstate of $\hat J_z=\hat L_z+\hat S_z$ with eigenvalue $M=m_\ell+m_s$. With this in mind organize your $\vert \ell m_\ell\rangle\vert sm_s\rangle$ states so that those with the same value of $M$ are on the same line. Explicitly, for instance, you would have
\begin{align}
\begin{array}{rlll}
M=\ell+s:&\vert \ell \ell\rangle\vert s s\rangle \\
M=\ell+s-1:& \vert \ell,\ell-1\rangle\vert ss\rangle&
\vert\ell\ell\rangle \vert s,s-1\rangle\\
M=\ell+s-2:&\vert \ell,\ell-2\rangle \vert ss\rangle & \vert\ell,\ell-1\rangle\vert s,s-1\rangle&\vert \ell \ell\rangle \vert s,s-2\rangle\\
\vdots\qquad& \qquad\vdots
\end{array}
\end{align}
and replace each state with a $\bullet$ to get
$$
\begin{array}{rlll}
\ell+s:&\bullet \\
\ell+s-1:&\bullet & \bullet\\
\ell+s-2:&\bullet&\bullet&\bullet\\
\vdots\qquad&\vdots
\end{array}
$$
Now, if $M=\ell+s$ is the largest value and it occurs once, the value of $j=\ell+s$ must occur once and also all the states $\vert j=\ell+s,m_j\rangle$ will occur once. There is a linear combination of the two states with $M=\ell+s-1$ that will be the state $\vert j=\ell+s,m_j=\ell+s-1\rangle$, there will be a linear combination of the three states with $M=\ell+s-2$ that will be the $\vert j=\ell+s,m_j=\ell+s-2\rangle$ state etc. Since we are only interested in enumerating the possible resulting values of $j$, and not interested in the actual states per se, we can eliminate from our table the first column since it contains one state with $m_j=\ell+s$, one with $m_j=\ell+s-1$ etc. Eliminating this column yields the reduced table
$$
\begin{array}{rll}
\ell+s-1: & \bullet\\
\ell+s-2:&\bullet&\bullet\\
\vdots\qquad&\vdots
\end{array}
$$
Since the value of $m_j=\ell+s-1$ occurs once, the value $j=\ell+s-1$ must occur once, and the states $\vert \ell+s-1,m_j\rangle$ will each occur once. We take out those from the list by deleting the first column to obtain a further reduced table
$$
\begin{array}{rl}
\ell+s-2:&\bullet\\
\vdots\qquad &\vdots
\end{array}
$$
The process so continues until exhaustion. In the examples above we have found $j=\ell+s,\ell+s-1$ and the final reduced table of the example, if not empty, would indicate the value of $j=\ell+s-1$. It is clear this process produces a decreasing sequence of $j$. The last value of $j$ is determined by the width of the original table. It is not hard to convince yourself that the width of the table will stop increasing once we reach $M=\vert \ell-s\vert $, and this is the last value of $j$. Thus by exhaustion you find the possible values of $j$ in the range
$$
\vert \ell-s\vert\le j\le \ell+s\, .
$$
As an example consider $\ell=1$ and $s=2$. The original table then looks like
$$
\begin{array}{rlll}
\frac{3}{2}:&\vert 11\rangle\vert 1/2,1/2\rangle \\
\frac{1}{2}:&\vert 10\rangle\vert 1/2,1/2\rangle & \vert 11\rangle\vert 1/2,-1/2\rangle\\
-\frac{1}{2}:&\vert 1,-1\rangle\vert 1/2,1/2\rangle&\vert 10\rangle\vert 1/2,-1/2\rangle\\
-\frac{3}{2}:&\vert 1,-1\rangle \vert 1/2,-1/2\rangle
\end{array}\qquad \to \qquad
\begin{array}{rlll}
\frac{3}{2}:&\bullet \\
\frac{1}{2}:&\bullet & \bullet \\
-\frac{1}{2}:&\bullet &\bullet \\
-\frac{3}{2}:&\bullet
\end{array}
$$
It is only $2$ column wide, and the width stop growing at $M=1/2$, indicating the possible $j$ in this case are $3/2$ and $1/2$, and indeed
$$
\vert 1-1/2\vert \le j\le 1+1/2
$$
Finally, note that the absolute value is required on the left because one could write state $\vert sm_s\rangle\vert \ell m_\ell\rangle$ without affecting the possible values of $j$.
It is true that the Dirac Hamiltonian does not commute with the spin vector $\vec{S}$. While we have $\left[\vec{J},H_{D}\right]=0$, the related commutator with $\vec{S}$ is nonzero. So the vector components—e.g. $S_{x}=\frac{\hbar}{2}\Sigma_{1}$ in terms of the $4\times4$ Dirac matrix $\Sigma_{1}$—do not represent conserved quantities.* However, the magnitude of the spin does commute with the Hamiltonian,
$$\left[\vec{S}\,^{2},H_{D}\right]=0.$$
So the magnitude of the spin is always $\frac{3}{4}\hbar^{2}$, and the spin is always described by a $s=\frac{1}{2}$ representation of the (universal cover of the) rotation group, $SU(2)$.
Note that an analogous commutator formula with $\vec{L}$ does not hold:
$$\left[\vec{L}\,^{2},H_{D}\right]\neq0.$$
This leads to the fact that the small and large components of a Dirac spinor describing an energy eigenvalue in a central potential have different eigenvalues of $l$, although they have the same value of $s$ and $j$.
*$H_{D}$ does commute with the helicity $\vec{S}\cdot\hat{p}$, but that is a nonlocal operator, since it involves a unit vector $\hat{p}=\vec{p}/|\vec{p}|$, where $\vec{p}=-i\hbar\vec{\nabla}$. (This is nonlocal because of the $1/|\vec{p}|$; functions of derivatives other than nonnegative integer powers cannot be defined locally.)
Best Answer
This problem appears interesting for the following reason. Let us write it down in Cartesian coordnates:
$$-\frac{1}{2}\left(\frac{\partial^2\psi}{\partial x^2}+\frac{\partial^2\psi}{\partial y^2}+\frac{\partial^2\psi}{\partial z^2}\right)+\frac{1}{2}(x\cdot S)^2\psi+L\cdot S\psi=E\psi$$
where I have introduced a 1/2 factor for later convenience. Now, I concentrate on x and I consider the operator
$$-\frac{1}{2}\frac{\partial^2}{\partial x^2}+\frac{1}{2}(x\cdot S)^2$$
One can introduce creation and annihilation operators in a similar way as for the harmonic oscillator
$$A_S=\frac{1}{\sqrt{2}}\left(\frac{\partial}{\partial x}+xS\right)$$
and the corresponding eigenvectors will be labeled as $|n,S\rangle$. The next step is to write down $L\cdot S=\frac{1}{2}(J^2-L^2-S^2)$ and we can restate this problem in the form
$$\left(A_S^\dagger A_S+\frac{1}{2}\right)\psi-\frac{1}{2}\left(\frac{\partial^2\psi}{\partial y^2}+\frac{\partial^2\psi}{\partial z^2}\right)+\frac{1}{2}(J^2-L^2-S^2)\psi=E\psi$$