Let us start from
$$
M_{ij}=\epsilon_{ijk}J_{k}\ ,\quad M_{0i}=-K_i
$$
(the minus sign is my convention, sorry!). If you know, a priori, that
\begin{align*}
[J_i,J_j]&=i\epsilon_{ijk}J_k\ ,\\
[K_i,K_j]&=-i\epsilon_{ijk}J_k\ ,\\
[J_i,K_j]&=i\epsilon_{ijk}K_j\ ,
\end{align*}
then you can just do a brute force computation for appropriate indices and figure out the pattern for the general case. For example
\begin{align*}
[M_{0i},M_{0j}]&=[K_i,K_j]=-i\epsilon_{ijk}J_k=-i\,M_{ij}\ ,
\end{align*}
and so on. Alternatively, you can write down $J$ explicitly in terms of $M$ by doing
$$
\epsilon_{ijq}M_{ij}=\epsilon_{ijq}\epsilon_{ijk}J_k=2\delta_{qk}J_k\Rightarrow J_k=\frac{1}{2}\epsilon_{kij}M_{ij}\ .
$$
Then we can substitute this directly on the commutation relation for $J$'s
$$
\frac{1}{4}[\epsilon_{mni}M_{mn},\epsilon_{pqj}M_{pq}]=i\epsilon_{ijk}\epsilon_{uvk}M_{uv}=i\left(\delta_{iu}\delta_{jv}-\delta_{iv}\delta_{ju}\right)M_{uv}=-2iM_{ij}\ .\quad (\Box)
$$
Then you must replace the product of Levi-Civita symbols on the LHS of $(\Box)$ by their representation of the product of $\delta$'s to get
$$
[M_{ik},M_{jk}]=-2iM_{ij}\ .
$$
You then carry on substituting $J$ on the commutation relation for $J$ with $K$, but in the end you still must guess the general pattern.
Unfortunately none of these approaches are very practical. If all you want is a way to find the commutation relations for the generators of the Lorentz group, I suggest writing them down on a particular representation, e.g.
$$
(M_{\mu\nu})^{\sigma}_{\ \ \rho}=i\left(\eta_{\mu\rho}\delta^\sigma_{\ \ \nu}-\eta_{\nu\rho}\delta^\sigma_{\ \ \mu}\right)\ .
$$
On the above expression, $\sigma$ and $\rho$ represent matrix indices sort of speak. Then one knows how $M_{\mu\nu}$ transforms under a Lorentz transformation
$$
\Lambda M_{\mu\nu}\Lambda^{-1}=M_{\lambda\sigma}\Lambda^{\lambda}_{\ \ \mu}\Lambda^{\sigma}_{\ \ \nu}.\quad (\star)
$$
Considering infinitesimal transformations $\Lambda=1-\frac{i}{2}\xi^{\lambda\sigma}M_{\lambda\sigma}$ (for some antisymmetric parameter $\xi$), the LHS of $(\star)$ can be written as
$$
\Lambda M_{\mu\nu}\Lambda^{-1}=M_{\mu\nu}+\frac{i}{2}\xi^{\lambda\sigma}[M_{\mu\nu},M_{\lambda\sigma}]+\mathcal{O}(\xi^2)\ ,
$$
while the RHS of $(\star)$ can be written as
$$
M_{\lambda\sigma}\Lambda^{\lambda}_{\ \mu}\Lambda^{\sigma}_{\ \nu}=M_{\mu\nu}-\frac{1}{2}\xi^{\lambda\sigma}\left(M_{\mu\lambda}\eta_{\nu\sigma}-M_{\mu\sigma}\eta_{\nu\lambda}+M_{\lambda\nu}\eta_{\mu\sigma}-M_{\sigma\nu}\eta_{\mu\lambda}\right)+\mathcal{O}(\xi^2)\ ,
$$
thus implying
$$
[M_{\mu\nu},M_{\lambda\sigma}]=i\left(M_{\mu\lambda}\eta_{\nu\sigma}-M_{\mu\sigma}\eta_{\nu\lambda}+M_{\lambda\nu}\eta_{\mu\sigma}-M_{\sigma\nu}\eta_{\mu\lambda}\right)\ .
$$
Hope this helps! Cheers.
