Quaternion matrices are matrices whose elements are quaternions, i.e. numbers of the form \(\alpha_ R+\alpha_ Ii+\beta_ Rj-\beta_ Ik\) where \(i^ 2=j^ 2=k^ 2=-1\), \(ij=-ji=k\), \(jk=-kj=i\), and \(ki=-ik=j\). Such matrices arise naturally in e.g. quantum mechanical problems. All the steps of the classical Francis QR-algorithm for computing the eigenvalues and vectors of a complex matrix have quaternion analogies. This paper describes all these steps and thus develops a quaternion QR- algorithm with implicit shifts which, by means of a sequence of quaternion unitary similarity transformations, produces a Schur-like triangular matrix. The diagonal elements of this matrix are representatives of the wanted eigenvalues of the matrix. Any quaternion \(n\times n\) matrix can be written in the form \(A+jB\) where both A and B are \(n\times n\) complex matrices. The algorithm works directly with the matrices A and B and preserves quaternion structure throughout. It is backward stable.