In this Wikipedia page it is said that the square roots of -1 in the quaternion ring are the elements of the imaginary sphere. I don't understand why this is so. I don't understand the system that's written there. If I do
$$(a+bi+cj+dk)^2=a^2-b^2-c^2-d^2+2(ab+cd)i+2(ac-bd)j+2(ad+bc)k=-1$$
the system we get is
$$\begin{cases} a^2-b^2-c^2-d^2=-1 \\ ab+cd=0 \\ ac-bd=0 \\ ad+bc=0 \end{cases}$$
I suppose this could be solved to get the same solution, but it is quite tiresome and I'd like to understand how to get to the prettier system the Wikipedia gives.
Add your first equation to $a^2+b^2+c^2+d^2=1$ and you get $a=0$, whence $b^2+c^2+d^2=1$.
Then observe
$$(b\mathbf{i}+c\mathbf{j}+d\mathbf{k})^2=-(b^2+c^2+d^2)+(cd-dc)\mathbf{i}+(db-bd)\mathbf{j}+(bc-cb)\mathbf{k} $$ $$=-1+0\mathbf{i}+0\mathbf{j}+0\mathbf{k}=-1.$$ Hence any quaternion $\mathbf{q}=a+b\mathbf{i}+c\mathbf{j}+d\mathbf{k}$ with $a=0,b^2+c^2+d^2=1$ is a square root of $-1$. The error in your calculations is that you didn't take into account anticommutativity, i.e. $$(c\mathbf{j})(d\mathbf{k})=cd\;\mathbf{i}\quad\text{but}\quad (d\mathbf{k})(c\mathbf{j})=-dc\;\mathbf{i}.$$