How to find invertible matrix $P$ Such that $B=P^{-1}AP$ 
Let $A,B$ be the two $3 \times 3$ matrices
$$A=\begin{bmatrix}
3 & 2 &-5\\2 &6&-10\\1 &2 & -3 \end{bmatrix}$$
$$B=
\begin{bmatrix}
6 &20 &-34\\6 &32&-51\\4 &20 &-32 \end{bmatrix}
$$
Suppose that there is a Non Singular matrix $P$ that $P^{-1}AP=B$. Find $P$

My idea is using Eigen Values we get Diagonilization of Matrices $A$ and $B$ as:
$$A=RD_AR^{-1}$$
$$B=QD_BQ^{-1}$$
Then we get:
$$QD_BQ^{-1}=P^{-1}RD_AR^{-1}P$$
But i feel this is very tedious procedure?
 A: Let $$P=\begin{bmatrix} a & b& c \\ u & v& w\\ x & y & z \end{bmatrix}$$
$$PB=AP.$$
Then set up 9 linear simultaneous  homogeneous equations for $a,b,...y,z$ as
$$C_{9 \times 9} X_{9\times 1} =O_{9\times 1}$$
many solutions will exist for $X$, from a vector $X$, one can then determine a matrix $P$. The matric $P$ will be non-unique. The matrix $C$ will have its $\det|C|=0$ as zero, otherwise $X=O.$
A: When we don't know explicitly the eigenvalues, there are two methods


*

*We solve the equation $PB=AP$; the space of solutions has dimension $dim(C(A))$ where $C(A)$ is the commutant of $A$ (here $5$). The general solution is in the form $P=a_1P_1+\cdots+a_5P_5$ where the $(P_i)$ are known matrices.  We randomly choose the coefficients $(a_i)$ and we obtain (with probability $1$) an invertible matrix $P$ and we are done (of course, if your matrix $P$ is not invertible, then you start again for free).

*We calculate the Frobenius normal forms (with a PC, it's very fast) $A=QFQ^{-1},B=RFR^{-1}$ (with the same $F$) and we deduce $P$.
Yet, here we know $spectrum(A)=\{2,2,2\}$ and that  $dim(\ker(A-2I))=2$. 
$e_1=(1,0,0)$ is not in $\ker(A-2I)$ or in $\ker(B-2I)$ and $e_2=(B-2I)(e_1)=[4,6,4]^T\in \ker(B-2I),e'_2=(A-2I)(e_1)=[1,2,1]^T\in\ker(A-2I)$.
We complete the bases of $\ker(B-2I),\ker(A-2I)$ with $e_3=[-5,1,0]^T$ and $e'_3=[-2,1,0]^T$.
We choose $P$ s.t. $P(e_1)=e_1,P(e_2)=e'_2,P(e_3)=e'_3$.
Finally $P\begin{pmatrix}1&4&-5\\0&6&1\\0&4&0\end{pmatrix}=\begin{pmatrix}1&1&-2\\0&2&1\\0&1&0\end{pmatrix}$ and $P=\begin{pmatrix}1&3&-21/4\\0&1&-1\\0&0&1/4\end{pmatrix}$.
A: If you want, you can do it by solving a single explicit 9 by 9 system (also tedious but straightforward!).
Indeed, call $p_i =Pe_i$. Note that $p_i^l =P_{il}$ cause we are using standard basis vectors. Here $p_i^l$ denotes the $l$ component of $p_i$. 
We have $A=PBP^{-1}$. Applying these matrices to $p_i$ and taking the k-th component we get:


*

*LHS: $(Ap_i)_k = \sum A_{kl} p_i^l$

*RHS: $(PBe_i)_k = (P\sum B_{ij} e_j)_k = (\sum B_{ij} Pe_j)_k = \sum B_{ij} p_j^k$
By equating LHS and RHS we get 9 equations (varying $i,k$) in 9 variables (the $p_i^l$'s). As far as I know, you should solve 6 linear system 3x3 to find the eigenvectors, so this should be faster.
A: This is an ad-hoc method that works for almost simple matrices. It is optimized for dimension 3. If you dont want to see the proof, skip at the end to see the recipe for computation.
Definition. A simple matrix is a matrix that can be written as $v w^T$ for some vectors $v,w$. In particular $(v w^T)_{ij} = v_i w_j$.
If you know what tensors are, I called it simple because it corresponds to a simple tensor $v \otimes w$ in the correspondence $Hom(W,V) \simeq W^* \otimes V$.
Observation. A matrix $M$ is simple if and only if the rows are all scalar multiples of the same vector. In this case, for a choice of $M_{11} = ab$, we have that $M =uv^T$ where $v$ is the first row divided by $a$, and $u$ is the first column divided by $b$.
First step: reduce to a simple matrix. Note that 
$$ A = 2I + uv^T, B= 2I + wz^T$$
where $u= (1\ \  2\ \  1), v = (1\ \ 2\ \ -5), w=(2 \ \  3\ 2), z=(2\ \ 10 \ \  -17) $. 
To check this computation, use the observation with $A_{11} = 1 \cdot 1, B_{11} = 2 \cdot 2$. Set $C= uv^T, D=wz^T$ for the corresponding simple matrices.
Suppose that we find a matrix $P$ such that $P^{-1}(uv^T)P = wz^T$. Then
$$ P^{-1}AP = P^{-1}( uv^T + 2I)P = P^{-1}(uv^T)P + 2P^{-1}P = wz^T+2I = B$$
So we would be done.
Second step: reduce to matrices with same frobenius norm.
Next paragraph will explain why this step is necessary. Recall that the frobenius norm $\| \cdot \|$ is the norm of the matrix seen as a vector, i.e. the square root of the sum of squared components.
Use a matrix of the form $X_{\lambda} = \text{diag}( \lambda\ \  1 \ \ 1)$ with $\lambda \neq 0$ to conjugate $C$. You will get 
$$\| X_{\lambda^{-1}} C X_{\lambda} \|^2 = C_{11}^2 + \lambda^2(C_{12}^2 + C_{13}^2) + \lambda^{-2}(C_{21}^2 + C_{31}^2) + (\sum_{i,j>1} C_{ij}^2) $$
Note that the frobenius norm of $D$ is strictly greater than the frobenius norm of $C$ (the expression above for $\lambda=1$). Moreover, this expression in lambda goes to infinity for $\lambda \to \infty$, so that there exist a $\lambda$ such that
$$\| X_{\lambda^{-1}}C X_{\lambda} \|^2 = \| D \|^2$$
You can find it by solving a quadratic equation.
Also, note that $X_{\lambda^{-1}}C X_{\lambda} = (X_{\lambda^{-1}}u)(X_{\lambda}v)^T$ is still a simple matrix. Set $\alpha = X_{\lambda^{-1}}u, \beta = X_{\lambda}v$
Third step: conjugating simple matrices with same frobenius norm.
Firstly, note that for a simple matrix
$$ \| u \|^2 \| v\|^2 = (\sum_i u_i^2 ) (\sum_j v_j^2 ) = \sum_{i,j} (u_iv_j)^2 = \| uv^T \|^2 $$
i.e. the frobenius norm is the product of the norms of the two vectors. 
We look for an orthogonal $Q$ such that $(Q\alpha)(Q^T\beta)^T = Q \alpha \beta^T Q^{-1} = wz^T$. It is enough if we find a $Q$ and a $\mu$ such that
$$ Q\alpha = \mu w, Q\beta = \mu^{-1} z$$
Taking norms and recalling that $Q$ is orthogonal:
$$ \| \alpha \| = \mu \| w\|, \| \beta \| = \mu^{-1} \| z \| $$
So that $\mu = \|\alpha\| / \|w\|$. Indeed, the second equation is satisfied by
$$\mu^{-1} \|z\| = \| w\| \| z\| / \|\alpha \| = \|wz^T \| / \|\alpha\| = \| \alpha \beta^T \| / \| \alpha \| = \| \beta\| $$
Now set $\gamma = \mu w, \delta = \mu^{-1} z$. Note that we have
$$\| \alpha\| = \|\gamma\|, \|\beta\| = \|\delta\|$$
$$ \langle \alpha, \beta \rangle = \text{tr}(\alpha \beta^T) = \text{tr}(uv^T) = \text{tr}(A) - \text{tr}(2I) = \text{tr}(B) -\text{tr}(2I) =\text{ tr}(wz^T) = \langle \gamma, \delta \rangle $$
There exist an orthogonal matrix $Q$ such that $Q \alpha = \gamma, Q \beta = \delta, Q( \alpha \times \beta) = Q(\gamma \times \delta)$. You can find this by an ordinary change of basis process. 
Recipe. On balance, the process is the following:


*

*As in first step, set $C=A-2I, D= B-2I$ and $u= (1\  2\  1), v = (1\ 2\ -5), w=(2 \ 3\ 2), z=(2\ 10 \ -17) $. 

*Solve the following biquadratic equation for $\lambda$:
$$ C_{11}^2 + \lambda^2(C_{12}^2 + C_{13}^2) + \lambda^{-2}(C_{21}^2+C_{31}^2) + \sum_{i,j >1} C_{ij}^2 = \sum_{i,j} D_{ij}^2 $$
And set $\alpha = (\lambda^{-1} \ 2 \ 1), \beta = (\lambda \ 2 \ -5)$.

*Set $\mu = \| \alpha \| / \| w\| $, and $\gamma = \mu w, \delta = \mu^{-1}z$.

*Find a base change $Q$ such that $Q \alpha = \gamma, Q \beta = \delta, Q (\alpha \times \beta) = \gamma \times \delta$.


Then your final matrix $P$ will be $P=Q^T \text{diag}(\lambda \ \ 1 \ \ 1)$.
Pros and cons: the pro is that you will just have to compute an explicit quadratic equation and a base change. The con is that probably big numbers will be involved!
