Let $\mathbb{H}$ be the real quaternion division ring consisting of all elements of form: $a+bi+cj+dk$ in which $a,b,c,d\in\mathbb{R},i^2=j^2=k^2=-1, ij=-ji=k$ and $\alpha=a+bi+cj+dk,\beta=a+\sqrt{b^2+c^2+d^2}i\in\mathbb{H}$.
Let us consider $x=a_1+b_1i+c_1j+d_1k\in\mathbb{H}$. Assume that $\alpha x=x\beta$. Then, $$(aa_1-bb_1-cc_1-dd_1)+(ab_1+ba_1+cd_1-dc_1)i+(ac_1+ca_1-bd_1+db_1)j+(ad_1+da_1+bc_1-cb_1)k$$
$$=(a_1a-b_1\sqrt{b^2+c^2+d^2})+(b_1a+a_1\sqrt{b^2+c^2+d^2})i+(c_1a+d_1\sqrt{b^2+c^2+d^2})j+(ad_1-c_1\sqrt{b^2+c^2+d^2})k.$$ Therefore, $$(aa_1-bb_1-cc_1-dd_1)=(a_1a-b_1\sqrt{b^2+c^2+d^2})$$ $$(ab_1+ba_1+cd_1-dc_1)=(b_1a+a_1\sqrt{b^2+c^2+d^2}),$$ $$(ac_1+ca_1-bd_1+db_1)=(c_1a+d_1\sqrt{b^2+c^2+d^2}),$$ $$(ad_1+da_1+bc_1-cb_1)=(ad_1-c_1\sqrt{b^2+c^2+d^2}),$$ that is, $$(-bb_1-cc_1-dd_1)=(-b_1\sqrt{b^2+c^2+d^2})$$ $$(ba_1+cd_1-dc_1)=(a_1\sqrt{b^2+c^2+d^2}),$$ $$(ca_1-bd_1+db_1)=(d_1\sqrt{b^2+c^2+d^2}),$$ $$(da_1+bc_1-cb_1)=(-c_1\sqrt{b^2+c^2+d^2}).$$ I would like to ask whether it is possible to continue.
Any counterexample or reference or technique is very much appreciated. Thank you in advance.
Note that what you get is an honest linear system, and you can just solve it. You'll get a $2$-dimensional space of solutions. I'll try to explain how to be smart about the problem (but you don't have to be, you can just solve your system).
You can start by cleaning things up to simplify computations. First, note that the real part of $\alpha$ and $\beta$ is the same and does not play any role: if we write $\alpha= a+\alpha_0$ and $\beta=a+\beta_0$ where $\alpha_0$ and $\beta_0$ are pure quaternions, then $\alpha x=x\beta$ is equivalent to $\alpha_0x=x\beta_0$, since scalars commute with everything. So we might as well assume that $a=0$ (you can check that indeed $a$ does not appear in your system).
Also, you can divide $\alpha$ and $\beta$ by their norm $N=\sqrt{b^2+c^2+d^2}$ so get $\alpha'=b'i+c'j+d'ij$ and $\beta'=i$. (If $N=0$, then $\alpha=\beta=0$ and all quaternions $x$ are solutions, so I will assume that this is not the case.) Obviously this changes nothing to the problem and yields the same solutions $x$, so in the end we can assume that $\alpha$ is a pure quaternion with unit norm, and $\beta=i$. (This is just for convenience, it does not really change anything.)
Now any quaternion $x$ has a unique decomposition $x=x_0+x_1$ where $x_0$ commutes with $i$ and $x_1$ anti-commutes with $i$ (explicitly, $x_0$ is in $\operatorname{Vect}(1,i)$ and $x_1$ is in $\operatorname{Vect}(j,ij)$). And it is easy to see that if $x$ and $x'$ are two solutions to your problem with $x\neq 0$ (of course $x=0$ is always a trivial solution), then there is a unique $q$ which commutes with $i$ such that $x'=xq$. So once you find some non-trivial solution $x$, the set of all solutions is $$ \{xq\,|\,q \text{ commutes with } i\} = \{x(\lambda +\mu i )\,|\, \lambda,\mu\in \mathbb{R} \}. $$
Now if $x=x_0+x_1$, then $x'=x_0q+x_1q$ and you can see that $x_0q$ commutes with $i$ and $x_1q$ anti-commutes with $i$, so $x_0'=x_0q$ and $x_1'=x_1q$. But since $\operatorname{Vect}(1,i)\simeq \mathbb{C}$ as an algebra, it means that you can choose $x_0'$ as you want, and then solve $x_0q=x_0'$ in $\mathbb{C}$, which of course has a unique solution $q=x_0'/x_0$ if $x_0\neq 0$ since $\mathbb{C}$ is a field.
So in the end there are two cases:
if $\alpha=-i$, then $x$ can be any quaternion which anti-commutes with $i$, so $x_0=0$ and you can choose $x_1$ to be anything in $\operatorname{Vect}(j,ij)$;
if $\alpha \neq -i$, then you can choose $x_0$ to by anything in $\operatorname{Vect}(1,i)\setminus \{0\}$, and then $x_1$ will be uniquely determined. For instance, you can choose $x_0=1$ if you want.
Now when $\alpha\neq -i$ you do have to find $x_1$, by solving the system you had in your question, but in a simpler version. To summarize, in terms of your system:
if $b=c=d=0$, then any $x$ is a solution;
if $c=d=0$ and $b<0$, then $a_1=b_1=0$ and $c_1$ and $d_1$ can be anything;
otherwise, you can solve your system by assuming that $a_1=1$ and $b_1=0$ (and replace $\sqrt{a^2+b^2+c^2}$ by $N$ to make things prettier), and then you will find a unique solution; all other solutions are of the form $x(\lambda+\mu i)$.