As a first step, I'd move $P_1$ to the origin, so that the points become $P_1=(0,0,0)$, $P_2=(x_2-x_1, y_2-y_1, z_2-z_1)$, $P_3=(x_3-x_1, y_3-y_1, z_3-z_1)$. From there it's just a question of applying a rotation matrix.
A rotation matrix $\mathbf{R}$ which preserves all distances and shapes etc is
given by
$$
\mathbf{R}=\left(
\begin{array}
[c]{ccc}%
e_{1}^{x} & e_{2}^{x} & e_{3}^{x}\\
e_{1}^{y} & e_{2}^{y} & e_{3}^{y}\\
e_{1}^{z} & e_{2}^{z} & e_{3}^{z}%
\end{array}
\right)
$$
where the vectors $\left( e_{1}^{\alpha},e_{2}^{\alpha},e_{3}^{\alpha
}\right) $ for $\alpha\in\left\{ x,y,z\right\} $ are orthogonal and of unit
length, i.e.
$$
\sum_{i=1}^{3}e_{i}^{\alpha}e_{i}^{\beta}=\left\{
\begin{array}
[c]{ccc}%
1 & & \alpha=\beta\\
0 & & \alpha\neq\beta
\end{array}
\right. .%
$$
Applying $\mathbf{R}$, your rotated points $P_{i}^{\prime}=\left(
x_{i}^{\prime},y_{i}^{\prime},z_{i}^{\prime}\right) $ are given by
$$
\left(
\begin{array}
[c]{c}%
x_{i}^{\prime}\\
y_{i}^{\prime}\\
z_{i}^{\prime}%
\end{array}
\right) =\mathbf{R}\left(
\begin{array}
[c]{c}%
x_{i}\\
y_{i}\\
z_{i}%
\end{array}
\right)
$$
I suggest making the top row of $\mathbf{R}$ the unit
vector in the direction from $P_{1}$ to $P_{2}$, i.e.
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{x}\\
e_{2}^{x}\\
e_{3}^{x}%
\end{array}
\right) =\frac{1}{\sqrt{\left( x_{1}-x_{2}\right) ^{2}+\left( y_{1}%
-y_{2}\right) ^{2}+\left( z_{1}-z_{2}\right) ^{2}}}\left(
\begin{array}
[c]{c}%
x_{1}-x_{2}\\
y_{1}-y_{2}\\
z_{1}-z_{2}%
\end{array}
\right)
$$
which would mean that the line from $P_{1}$ to $P_{2}$ gets mapped into the
$x$-axis.
Working out what you could use for the second and third rows of $\mathbf{R}$ is now just a question of solving some simple linear equations, to make sure that the line from $P_{1}$ to $P_{3}$ has no z component.
To solve for $\left( e_{1}^{z},e_{2}^{z},e_{3}^{z}\right)$, the vector
$e_{i}^{z}$ must have a zero dot-product with both $e_{i}^{x}$ and $\left(
x_{3}-x_{1},y_{3}-y_{1},z_{3}-z_{1}\right) $, so the equations are
$$
\begin{align*}
e_{1}^{x}e_{1}^{z}+e_{2}^{x}e_{2}^{z}+e_{3}^{x}e_{3}^{z} & =0\\
\left( x_{3}-x_{1}\right) e_{1}^{z}+\left( y_{3}-y_{1}\right) e_{2}%
^{z}+\left( z_{3}-z_{1}\right) e_{3}^{z} & =0
\end{align*}
$$
Hence
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{z}\\
e_{2}^{z}\\
e_{3}^{z}%
\end{array}
\right) =\lambda_{z}\left(
\begin{array}
[c]{c}%
\left( z_{3}-z_{1}\right) e_{2}^{x}-\left( y_{3}-y_{1}\right) e_{3}^{x}\\
\left( x_{3}-x_{1}\right) e_{3}^{x}-\left( z_{3}-z_{1}\right) e_{1}^{x}\\
\left( y_{3}-y_{1}\right) e_{1}^{x}-\left( x_{3}-x_{1}\right) e_{2}^{x}%
\end{array}
\right)
$$
for some $\lambda_z$, which should be determined so that $e_{i}^{z}$ is a vector
of unit length. Finally $e_{i}^{y}$ can be determined as the vector which is
perpendicular to both $e_{i}^{x}$ and $e_{i}^{z}$, i.e.
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{y}\\
e_{2}^{y}\\
e_{3}^{y}%
\end{array}
\right) =\lambda_{y}\left(
\begin{array}
[c]{c}%
e_{3}^{z}e_{2}^{x}-e_{2}^{z}e_{3}^{x}\\
e_{1}^{z}e_{3}^{x}-e_{3}^{z}e_{1}^{x}\\
e_{2}^{z}e_{1}^{x}-e_{1}^{z}e_{2}^{x}%
\end{array}
\right)
$$
where again $\lambda_{y}$ is determined so that $e_{i}^{y}$ is a vector of
unit length.
Best Answer
A distance-preserving map is known as an isometry. Fix $x_0 \in \Pi_1$ and $y_0 \in \Pi_2$. To find an isometry between $\Pi_1$ and $\Pi_2$, mapping $x_0$ to $y_0$, you can follow the following steps:
In fact, this method will generate all possible isometries between these planes, as the Mazur-Ulam theorem implies.
EDIT: $T$ is a linear transformation between two-dimensional spaces $\Pi'_1$ and $\Pi'_2$. Given the two orthonormal bases for $\Pi'_1$ and $\Pi'_2$ found in step 2, the matrix for the transformation, from one orthonormal basis to the other, is the identity matrix. That is, the transformation maps a coordinate vector in one basis to exactly the same coordinate vector in the other basis.
It is possible (and probably easier computationally) to extend $T$ to a linear isometry $S : \mathbb{R}^3 \to \mathbb{R}^3$, which will have a standard matrix.
To do this, extend the bases for $\Pi'_1$ and $\Pi'_2$ to bases for $\mathbb{R}^3$, then apply Gram-Schmidt. The result should be two orthonormal bases for $\mathbb{R}^3$. Again, define a linear map that maps one basis to the other. The result is a linear isometry that is the same as $T$ when restricted to $\Pi'_1$.
Let's go through an example. Suppose $\Pi_1$ is the plane defined by the equation $x + 2z + 3 = 0$ and $\Pi_2$ is the plane defined by $-x + 3y + z - 4 = 0$. In particular, fix points $x_0 = (-1, -1, -1) \in \Pi_1$ and $y_0 = (4, 0, 0) \in \Pi_2$.
Step 1: \begin{align*} \Pi_1' &= \lbrace (x, y, z) \in \mathbb{R}^3 : x + 2z = 0 \rbrace \\ \Pi_2' &= \lbrace (x, y, z) \in \mathbb{R}^3 : -x + 3y + z = 0 \rbrace \\ \end{align*}
Step 2: We find arbitrary, non-parallel vectors in $\Pi_1'$ and $\Pi_2'$. I'm going to take $(-2, 0, 1), (0, 1, 0) \in \Pi_1'$ and $(3, 1, 0), (1, 0, 1) \in \Pi_2'$. Since we also want to find $S$, we need one more vector each, that doesn't lie in the respective planes. For both, I'm going to choose $(1, 0, 0) \notin \Pi_1' \cup \Pi_2'$. We therefore start with (not orthonormal) bases: \begin{align*} \mathcal{B}_1 &= ((-2, 0, 1), (0, 1, 0), (1, 0, 0)) \\ \mathcal{B}_2 &= ((3, 1, 0), (1, 0, 1), (1, 0, 0)). \end{align*} We need to apply the Gram-Schmidt procedure to each of these bases. For $\mathcal{B}_1$, we have \begin{align*} u_1 &= \frac{(-2, 0, 1)}{\|(-2, 0, 1)\|} = \left(\frac{-2}{\sqrt{5}}, 0, \frac{1}{\sqrt{5}}\right) \\ u_2 &= \frac{(0, 1, 0) - ((0, 1, 0) \cdot u_1)u_1}{\|(0, 1, 0) - ((0, 1, 0) \cdot u_1)u_1\|} = (0, 1, 0) \\ u_3 &= \frac{(1, 0, 0) - ((1, 0, 0) \cdot u_1)u_1 - ((1, 0, 0) \cdot u_2)u_2}{\|(1, 0, 0) - ((1, 0, 0) \cdot u_1)u_1 - ((1, 0, 0) \cdot u_2)u_2\|} = \left(\frac{1}{\sqrt{5}}, 0, \frac{2}{\sqrt{5}}\right) \end{align*} So, we get an orthonormal basis $$\mathcal{E}_1 = \left(\left(\frac{-2}{\sqrt{5}}, 0, \frac{1}{\sqrt{5}}\right), (0, 1, 0), \left(\frac{1}{\sqrt{5}}, 0, \frac{2}{\sqrt{5}}\right)\right).$$ Similarly, using $\mathcal{B}_2$, we obtain $$\mathcal{E}_2 = \left(\left(\frac{3}{\sqrt{10}}, \frac{1}{\sqrt{10}}, 0\right), \left(\frac{1}{\sqrt{110}}, \frac{-3}{\sqrt{110}}, \frac{10}{\sqrt{110}}\right), \left(\frac{1}{\sqrt{11}}, \frac{-3}{\sqrt{11}}, \frac{-1}{\sqrt{11}}\right)\right).$$ If we remove the last vector in $\mathcal{E}_1$ and $\mathcal{E}_2$, then we obtain orthonormal bases for $\Pi'_1$ and $\Pi'_2$ respectively.
Step 3: We define $S$ to be a linear transformation that maps \begin{align*} \left(\frac{-2}{\sqrt{5}}, 0, \frac{1}{\sqrt{5}}\right) &\mapsto \left(\frac{3}{\sqrt{10}}, \frac{1}{\sqrt{10}}, 0\right) \\ (0, 1, 0) &\mapsto \left(\frac{1}{\sqrt{110}}, \frac{-3}{\sqrt{110}}, \frac{10}{\sqrt{110}}\right) \\ \left(\frac{1}{\sqrt{5}}, 0, \frac{2}{\sqrt{5}}\right) &\mapsto \left(\frac{1}{\sqrt{11}}, \frac{-3}{\sqrt{11}}, \frac{-1}{\sqrt{11}}\right). \end{align*} To compute the standard matrix, we compute the image of the standard basis vectors. For example, we have \begin{align*}(1, 0, 0) &= \left((1, 0, 0) \cdot \left(\frac{-2}{\sqrt{5}}, 0, \frac{1}{\sqrt{5}}\right)\right)\left(\frac{-2}{\sqrt{5}}, 0, \frac{1}{\sqrt{5}}\right) \\ &+ ((1, 0, 0) \cdot (0, 1, 0))(0, 1, 0) \\ &+ \left((1, 0, 0) \cdot \left(\frac{1}{\sqrt{5}}, 0, \frac{2}{\sqrt{5}}\right)\right)\left(\frac{1}{\sqrt{5}}, 0, \frac{2}{\sqrt{5}}\right) \\ &= \frac{-2}{\sqrt{5}}\left(\frac{-2}{\sqrt{5}}, 0, \frac{1}{\sqrt{5}}\right) + 0(0, 1, 0) + \frac{1}{\sqrt{5}}\left(\frac{1}{\sqrt{5}}, 0, \frac{2}{\sqrt{5}}\right) \\ &\mapsto \frac{-2}{\sqrt{5}}\left(\frac{3}{\sqrt{10}}, \frac{1}{\sqrt{10}}, 0\right) + 0\left(\frac{1}{\sqrt{110}}, \frac{-3}{\sqrt{110}}, \frac{10}{\sqrt{110}}\right) + \frac{1}{\sqrt{5}}\left(\frac{1}{\sqrt{11}}, \frac{-3}{\sqrt{11}}, \frac{-1}{\sqrt{11}}\right) \\ &= \left(\frac{-33\sqrt{2} + \sqrt{55}}{55}, \frac{-11\sqrt{2} - 3\sqrt{55}}{55}, \frac{-\sqrt{55}}{55}\right). \end{align*} We have already $$(0, 1, 0) \mapsto \left(\frac{1}{\sqrt{110}}, \frac{-3}{\sqrt{110}}, \frac{10}{\sqrt{110}}\right)$$ Using similar methods, we finally have, $$(0, 0, 1) \mapsto \left(\frac{33\sqrt{2}-4\sqrt{55}}{110}, \frac{11\sqrt{2}-6\sqrt{55}}{110}, \frac{4\sqrt{55}}{110}\right).$$ This gives us the standard matrix for $S$: $$\mathcal{M}(S) = \begin{pmatrix} \frac{-33\sqrt{2} + \sqrt{55}}{55} & \frac{1}{\sqrt{110}} & \frac{33\sqrt{2} - 4\sqrt{55}}{110} \\ \frac{-11\sqrt{2} - 3\sqrt{55}}{55} & \frac{-3}{\sqrt{110}} & \frac{-11\sqrt{2} - 6\sqrt{55}}{110} \\ \frac{-\sqrt{55}}{55} & \frac{10}{\sqrt{110}} & \frac{4\sqrt{55}}{110} \end{pmatrix}.$$
Step 4: We define the isometry, mapping $x_0$ to $y_0$ by $$\begin{pmatrix} x \\ y \\ z \end{pmatrix} = \begin{pmatrix} \frac{-33\sqrt{2} + \sqrt{55}}{55} & \frac{1}{\sqrt{110}} & \frac{33\sqrt{2} - 4\sqrt{55}}{110} \\ \frac{-11\sqrt{2} - 3\sqrt{55}}{55} & \frac{-3}{\sqrt{110}} & \frac{-11\sqrt{2} - 6\sqrt{55}}{110} \\ \frac{-\sqrt{55}}{55} & \frac{10}{\sqrt{110}} & \frac{4\sqrt{55}}{110} \end{pmatrix}\left(\begin{pmatrix} x \\ y \\ z \end{pmatrix} - \begin{pmatrix} -1 \\ -1 \\ -1 \end{pmatrix}\right) + \begin{pmatrix} 4 \\ 0 \\ 0 \end{pmatrix}.$$