If one were to list 5 computational problems that are central to computer science, solving the matrix equation $Ax = b$ would be one of them. Here, $A$ is your matrix, $b$ is your result, and $x$ is the vector you are trying to solve for.
Regardless of what method you decide to use, you should be using a library and not manually coding the algorithm. There are two classical methods for solving the linear system $Ax = b$. The first is Gaussian elimination, suitable for small to medium-sized, dense matrices. This is implemented in NumPy as numpy.linalg.solve, or more generally, in LAPACK as gesv. Gaussian elimination, in all its implementations, are exact in exact arithmetic, and so you can never beat the $O(n^3)$ scaling.
The second is the iterative method, suitable for large, sparse matrices. There are several variations on this method. If $A$ is symmetric and positive-definite, you can use Conjugate Gradients. Otherwise, you will want to use GMRES. Actual library implementations of these are slightly more varied, and you will have to search the internet to find one suitable for your purposes. Iterative methods take advantage of sparsity to perform far faster than Gaussian elimination, at the cost of not being exact in exact arithmetic. However, given well-conditioned and sparse matrices, iterative methods have superb accuracy.
There is a third but less common option. If you know some special knowledge about the structure of your matrix, then you can also look for structured matrix solvers. Given the nature of your question, it is likely that this is not the case for your problem.
The bottom line is that you should search the internet to find a suitable library. If computational efficiency is very important to you, nothing that you come up with is going to beat the performance of subroutines in libraries like LAPACK.
Given a quaternion $q$ and a complex number $\lambda$, for scalar multiplication you can either apply $\lambda$ to $q$ from the left (i.e. $\lambda q$) or from the right (i.e. $q\lambda$). This means you can interpret $\mathbb{H}$ as a left complex vector space or as a right complex vector space. Now consider the left and right multiplication maps
$$ L_p(x)=px, \qquad R_p(x)=xp. \tag{1}$$
Then $L_p(x\lambda)=L_p(x)\lambda$ for all complex numbers, so $L_p$ is a linear transformation of $\mathbb{H}$ as a right complex vector space. And $R_p(\lambda x)=\lambda R_p(x)$, so $R_p$ is a linear transformation of $\mathbb{H}$ as a left complex vector space.
However, $L_p$ is not linear if we treat $\mathbb{H}$ as a left complex vector space, and $R_p$ is not linear if we treat $\mathbb{H}$ as a right vector space. This is because $\mathbb{H}$ is not commutative. (Exercise.)
Using $\{1,\mathbf{j}\}$ as a basis for $\mathbb{H}$ as a complex vector space, you would write an arbitrary quaternion as $z+w\mathbf{j}$ if you're thinking left vector space, and write $z+\mathbf{j}w$ if you're thinking right vector space. (Also notice $w$ comes before $z$ in the alphabet, so we're backwards alphabetically.) You're thinking left vector space but you're examining the right-linear transformation $L_p$. That's problematic.
If you look at $L_p$ and think right vector space, you should get
$$ \begin{bmatrix} z & -\overline{w} \\ w & \phantom{-}\overline{z} \end{bmatrix}. \tag{2} $$
(Exercise.)
Representing quaternions as matrices using $R_p$ is problematic since $R_p\circ R_q\ne R_{pq}$. Indeed, this would be a function $\mathbb{H}\to M_2(\mathbb{C})$ which is not a homomorphism but rather an anti-homomorphism, i.e. it satisfies $R_p\circ R_q=R_{qp}$ (the order of multiplication is reversed). In order to turn it into a homomorphism proper, one must either (pre)compose it with an anti-automorphism of $\mathbb{H}$ or (post)compose with an anti-automorphism of $M_2(\mathbb{C})$. Quaternion conjugation satisfies $\overline{pq}=\overline{q}\,\overline{p}$ and matrix transpose satisfies $(AB)^T=B^TA^T$, so we can use these as anti-automorphisms.
In the case of quaternion conjugation, we can take $p=z+w\mathbf{j}$, then its conjugate $\overline{p}=\overline{z}-w\mathbf{j}$, then the matrix of $R_{\overline{p}}$ you can calculate (exercise) to be
$$ \begin{bmatrix} \overline{z} & \overline{w} \\ -w & z \end{bmatrix}, \tag{3} $$
and if instead you just calculated $R_p$'s matrix and took the transpose you'd get $(2)$ again.
Complex conjugation (applied entry-wise to a matrix) is an automorphism of $M_2(\mathbb{C})$, i.e. it satisfies $\overline{AB}=\overline{A}\,\overline{B}$, so we may compose it with any of the representations above to get another valid representation.
Best Answer
Check your computation of $P_2$.
When computing $P_2$, instead of adding $A$ to $B$ to get $A+B$, you multiple $A$ and $B$ and get $A*B$ which is incorrect.