As the other two replies mentioned, the answer is no, spin-orbit coupling cannot be derived from the Pauli equation. Why? In a sense you could just say the term has been neglected in an approximation already made to arrive at the Pauli equation. The starting point you need is the Dirac equation.
To my eyes, the most convenient way to find the spin-orbit term is to perform the Foldy-Wouthuysen (FW) transformation on the Dirac equation including a central Coulomb potential. My favorite reference for this, albeit not a very commonly used one, is Chapter 3.4-3.6 of Brown & Carrington, Rotational Spectroscopy of Diatomic Molecules. I'll just roughly reproduce their derivation here, but I've heard that a pdf of that book is not hard to find...
Start with the Dirac Hamiltonian (in old-fashioned notation)
$${\mathcal H}= \beta m c^2 - e \phi + c \vec{\alpha} \cdot \vec{\pi}$$
$$ \equiv \beta m c^2 + V + O, $$
where $\vec{\pi}=\vec{P} + e \vec{A}$ and the operators $V=-e \phi$ and $O=\vec{\alpha}\cdot \vec{\pi}$ are even and odd with respect to the "big" and "small" parts of the electron wavefunction. The general idea of the Foldy-Wouthuysen transformation is to choose a unitary operator $e^{i {\mathcal S}}$ such that when we apply a unitary transformation
$$ {\mathcal H}' = e^{i {\mathcal S}} {\mathcal H} e^{-i {\mathcal S}} + i \hbar (\partial_t \, e^{i {\mathcal S}} ) e^{-i {\mathcal S}},$$
the odd parts of the transformed Hamiltonian are suppressed by powers of $c$. It turns out that this is possible, and the right choice of ${\mathcal S}$ is
$$ {\mathcal S} = - \frac{i \beta O}{2 m c^2}. $$
It so happens that you actually have to apply the FW transformation twice in order to suppress $O$ by two powers of $c$. After doing this, one of the resulting terms in the twice-transformed Hamiltonian is
$$ \frac{1}{8 m^2 c^2} \left[ \vec{\alpha} \cdot \vec{\pi} , \left[ \vec{\alpha} \cdot \vec{\pi}, (-e \phi) \right] + i \hbar \partial_t (\vec{\alpha} \cdot \vec{\pi}) \right] $$
$$ = -\frac{e \hbar^2}{8 m^2 c^2} \nabla \cdot \vec{\mathcal E} - \frac{g_S \mu_B}{2 m c^2} \vec{S} \cdot (\vec{\mathcal E} \times \vec{\pi}), $$
where $\vec{\mathcal E}=-\nabla \phi$ is the Coulomb electric field. The electric field $\vec{\mathcal E}$ can be replaced by the explicit Coulomb expression $\sim Z e^2 \vec{r}/r^3$, and the numerator $\vec{r} \times \vec{p}$ is precisely the orbital angular momentum $\vec{L}$. Thus we have a term in the Hamiltonian proportional to $\vec{S} \cdot \vec{L}$, which is spin-orbit coupling. Oh and by the way, the first term above I believe is the Darwin term proportional to $\delta^3(\vec{r})$.
The Foldy-Wouthuysen transformation can also be used to derive the Pauli term $\sim \vec{\sigma} \cdot \vec{\mathcal B}$ and the familiar leading order relativistic correction to the kinetic energy, proportional to $p^4$. For a free particle, it can also be shown to account for the phenomenon of zitterbewegung.