Square roots of $-1$ in quaternion ring

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}.$$

They have it right. You have not been careful enough in $\pm$ signs. The coefficient of $i$ does begin with $2ab$ as you say, but then we get $$ c j d k + d k c j = c d i - c d i=0.$$ And so forth.