Let $F$ be the center of the division ring $D$. Then $D$ represents its class in the Brauer group $Br(F)$. The opposite ring $D^{opp}$ represents the inverse element. The reason why for example the quaternions are isomorphic to their opposite algebra is that the quaternions are an element of order 2 in $Br(\mathbf{R})$, and hence equal to its own inverse in the Brauer group.
To get a division algebra that is not isomorphic to its opposite algebra we can use an element of order 3 in the Brauer group. One method for constructing those is to start with a (cyclic) Galois extension of number fields $E/F$ such that $[E:F]=3$. Let $\sigma\in Gal(E/F)$ be the generator. Let $\gamma\in F$ be an element that cannot be written in the form $\gamma=N(x)$, where $N:E\rightarrow F, x\mapsto x\sigma(x)\sigma^2(x)$ is the relative norm map. Consider the set of matrices
$$
\mathcal{A}(E,F,\sigma,\gamma)=\left\{
\left(\begin{array}{rrr}
x_0&\sigma(x_2)&\sigma^2(x_1)\\
\gamma x_1&\sigma(x_0)&\sigma^2(x_2)\\
\gamma x_2&\gamma\sigma(x_1)&\sigma^2(x_0)
\end{array}\right)\mid x_0,x_1,x_2\in E\right\}.
$$
A theorem of A. Albert tells us that this forms a division algebra with center $F$, and its order in $Br(F)$ is 3, so it will not be isomorphic to its opposite algebra. The theory is described for example in ch. 8 of Jacobson's Basic Algebra II. The buzzword 'cyclic division algebra' should give you some hits.
For a concrete example consider the following. Let $F=\mathbf{Q}(\sqrt{-3})$ and let
$E=F(\zeta_9)$, with $\zeta_9=e^{2\pi i/9}$. Then $E/F$ is a cubic extension of cyclotomic fields, $\sigma:\zeta_9\mapsto\zeta_9^4$. I claim that the element $2$ does not belong the image of the norm map. This follows from the fact 2 is totally inert in the extension tower $E/F/\mathbf{Q}$. Basically because $GF(2^6)$ is the smallest finite field of characteristic 2 that contains a primitive ninth root of unity. Now, if $2=N(x)$ for some $x\in E$, then 2 must appear as a factor (with a positive coefficient) in the fractional ideal generated by $x$. But the norm map then multiplies that coefficient by 3, and as there were no other primes above 2, we cannot cancel that. Sorry, if this is too sketchy.
Anyway (see Jacobson again), the product of the $\gamma$ elements modulo $N(E^*)$ is the operation in the Brauer group $Br(E/F)\le Br(F).$ Therefore the opposite algebra should correspond to the choice $\gamma=1/2$, (or to the choice $\gamma=4$, as $2\cdot4=N(2)$). So
$$
\mathcal{A}(E,F,\sigma,2)^{opp}\cong
\mathcal{A}(E,F,\sigma,1/2).
$$
and the choices $\gamma=2$ and $\gamma=1/2$ yield non-isomorphic division algebras, as their ratio $=4$ is not in the image of the norm map.
======================================================
Edit (added more details here, because another answer links to this answer): In general the cyclic division algebra construction works much the same for any cyclic extension $E/F$. When $[E:F]=n$ we get a set of $n\times n$ matrices with entries in $E$. The number $\gamma$ appears in the lower diagonal part of the matrix. The condition for this to be a division algebra is that $\gamma^k$ should not be a norm for any integer $k, 0<k<n$. Obviously it suffices to check this for (maximal) proper divisors of $n$. In particular, if $n$ is a prime, then it suffices to check that $\gamma$ itself is not a norm.
Edit^2: Matt E's answer here linear algebra over a division ring vs. over a field gives a simpler cyclic division algebra of $3\times3$ matrices with entries in the real subfield of the seventh cyclotomic field.
One is tempted (as I was originally), to argue as follows: since $M\cong M\oplus M\cong M\times M$, we have
$$R = \mathrm{Hom}(M,M) \cong \mathrm{Hom}(M,M\times M) \cong \mathrm{Hom}(M,M)\times\mathrm{Hom}(M,M) = R\times R,$$
and stop there. The problem with the argument is that in some instances we are dealing with these objects as $\mathbb{Z}$-modules rather than as $R$-modules (specifically the isomorphisms are isomorphism from $\mathrm{Hom}(M,M)$ to $\mathrm{Hom}(M,M\times M)$ is, right now, just an isomorphism as $\mathbb{Z}$-modules), so some care needs to be exercised to make the argument actually work.
First, note that $M\cong M\times M$ as $\mathbb{Z}$-modules.
Define a homomorphism $\varphi\colon M\to M\times M$ by maping $(m_i)$ to $\bigl((m_{2i-1}), (m_{2i})\bigr)$. That is, $(m_1,m_2,m_3,\ldots)$ maps to $\bigl( (m_1,m_3,m_5,\ldots),(m_2,m_4,m_6,\ldots)\bigr)$.
The elements of $R$ can be thought of as infinite "column-finite matrices"; that is, each endomorphism $M\to M$ corresponds to a family of functions $(\mathbf{f}_i)_{i\in\omega}$, with $\mathbf{f}_i\colon\mathbb{Z}\to M$; hence $\mathbf{f}_i = \sum_{j\in\omega}f_{ji}e_j$, where $e_j$ is the element of $M$ that has $1$ in the $j$th coordinate and $0$s elsewhere, $f_{ji}\in\mathbb{Z}$, and $f_{ji}=0$ for almost all $j$.
If we take an element of $\mathbf{f}=(f_{ij})$ of $R$, and compose it with the isomorphism $M\to M\times M$, we obtain a homomorphism $M\to M\times M$, given by $\bigl( (f_{i,2j-1}), (f_{i,2j})\bigr)$. So the map $\mathrm{Hom}(M,M)\to\mathrm{Hom}(M,M)\times \mathrm{Hom}(M,M)$ maps $(f_{ij})$ to $\bigl( (f_{i,2j-1}),(f_{i,2j})\bigr)$.
If we compose maps $\mathbf{f}=(f_{ij})$ and $\mathbf{g}=(g_{ij})$, we get the map
$(h_{ij})$, where
$$h_{ij} = \sum_{k=1}^{\infty}g_{ik}f_{kj}.$$
Since $f_{kj}=0$ for almost all $k$, the sum is finite and $h_{ij}$ makes sense.
This gives the action of $R$ on $\mathrm{Hom}(M,M)$. The action of $R$ on $\mathrm{Hom}(M,M)\times\mathrm{Hom}(M,M)$ is then given by
$$
\mathbf{g}\bigl((f_{ij}),(f'_{ij})\bigr) = \bigl( \mathbf{g}(f_{ij}), \mathbf{g}(f'_{ij})\bigr)
= \bigl( (h_{ij}), (h'_{ij})\bigr)$$
where
$$h_{ij} = \sum_{k=1}^{\infty}g_{ik}f_{kj},\quad\text{and}\quad h'_{ij}=\sum_{k=1}^{\infty}g_{ik}f'_{kj}.$$
To see that the map from $\mathrm{Hom}(M,M)$ to $\mathrm{Hom}(M,M\times M)$ respects the action of $R$ (that is, that we get an $R$-module homomorphism, not merely a $\mathbb{Z}$-module homomorphism), let $(f_{ij})\in\mathrm{Hom}(M,M)$ and $\mathbf{g}\in R$. The element $\mathbf{g}(f_{ij})$ is mapped to
$$\left( \Bigl(\sum_{k=1}^{\infty}g_{ik}f_{k,2j-1}\Bigr),\Bigl( \sum_{k=1}^{\infty}g_{ik}f_{k,2j}\Bigr)\right)$$
On the other hand, if we let $\mathbf{g}$ act on $\bigl( (f_{i,2j-1}), (f_{i,2j})\bigr)$, we get
$$\mathbf{g}\left( (f_{i,2j-1}), (f_{i,2j})\right) = \left(\Bigl( \sum_{k=1}^{\infty}g_{ik}f_{k,2j-1}\Bigr), \Bigl( \sum_{k=1}^{\infty}g_{ik}f_{k,2j}\Bigr)\right),$$
that is, the same as the image of $\mathbf{g}(f_{ij})$. So the homomorphism $\mathrm{Hom}(M,M)\to\mathrm{Hom}(M,M)\times\mathrm{Hom}(M,M)$ given by $(f_{ij})\longmapsto ((f_{i,2k-1}),(f_{i,2k}))$ is actually a homomorphism as $R$-modules. Since the original map as $\mathbb{Z}$-modules was a bijection, so is this one, hence $\mathrm{Hom}(M,M)$ and $\mathrm{Hom}(M,M)\times\mathrm{Hom}(M,M)$ are isomorphic not only as $\mathbb{Z}$-modules, but also as $R$-modules. So we obtain:
$$R = \mathrm{Hom}(M,M) \cong \mathrm{Hom}(M,M)\times\mathrm{Hom}(M,M) =R\times R,$$
with the isomorphism being an isomorphism of $R$-modules, as desired.
Best Answer
The endomorphism ring of $R$ as a left $R$-module (equivalently, of $R^{op}$ as a right $R^{op}$-module) is $R^{op}$, which need not be isomorphic to $R$.
The correct statement for $R$ is then that $R$ is naturally isomorphic to the endomorphism of $R$ as a right $R$-module. This is Cayley's theorem for rings, and is a corollary of the (enriched) Yoneda lemma (although it also admits a straightforward direct proof).