Let $u$ and $v$ be vectors. Suppose the matrix $M$ is an isometry (preserves distances), takes $u$ to the $y$-axis and $v$ into the $x{-}y$ plane. Then it must take $u\times v$ to the $z$-axis and $u\times(u\times v)$ to the $x$-axis. That is
$$
M\begin{bmatrix}u&u\times v&u\times(u\times v)\end{bmatrix}=\begin{bmatrix}0&0&|u||u\times v|\\|u|&0&0\\0&|u\times v|&0\end{bmatrix}\tag{1}
$$
So we get
$$
M=\begin{bmatrix}0&0&|u||u\times v|\\|u|&0&0\\0&|u\times v|&0\end{bmatrix}\begin{bmatrix}u&u\times v&u\times(u\times v)\end{bmatrix}^{-1}\tag{2}
$$
Thus, $M$ takes three orthogonal vectors to three orthogonal vectors, preserving their lengths, so $M$ is an isometry.
If we let $u=b-a$ and $v=c-a$ and compute $M$ as in $(2)$, then $M(b-a)$ lies on the $y$-axis, and $M(c-a)$ lies in the $x{-}y$ plane. The only thing we can't do is map $a$ to the origin. A linear map sends $0$ to $0$, and since $M$ preserves distances $|Ma|=|M(a-0)|=|a-0|=|a|$. So unless $a=0$, $Ma\not=0$. To map $a$ to the origin, we need to apply a translation as well.
Compute $M$ as above, then the affine-linear map $A(x)=M(x-a)$ satisfies the requirements.
I would guess that you're supposed to use "homogeneous coordinates", and that the "3D" matrix you're supposed to compute is the matrix of a projective transformation. When we use homogeneous coordinates, a coordinate vector for a point in a plane has three components rather than two. This may seem unintuitive or unappealing, but there are advantages to homogeneous coordinates: 1) "points at infinity" can be assigned coordinates as easily as ordinary points ; 2) every projective transformation is given by multiplication by a $3 \times 3$ matrix. (There are other advantages too; mathematicians have discovered that homogeneous coordinates are a natural and elegant way to assign coordinates to points in a projective plane.)
Suppose we have a projective transformation that is represented by a matrix $H \in \mathbb R^{3 \times 3}$, so that point with coordinate vector $x$ maps to a point with coordinate vector $Hx$. If we have four points $x_1,\ldots,x_4$ in "general position" and we know these four points map to $y_1,\ldots, y_4$, then we can use this information to compute $H$.
We cannot simply say that $y_i = H x_i \, \,$, because one quirky thing about homogeneous coordinates is that if $y$ is a coordinate vector for a point in a projective plane, then any nonzero scalar multiple of $y$ is also a coordinate vector for the same point. This might seem weird but it's something we must get used to.
The most we can say is that $y_i$ is a (nonzero) scalar multiple of $H x_i$. An equivalent way to say the same thing is
\begin{equation}
y_i \times H x_i = 0
\end{equation}
where $\times$ denotes the cross product.
Notice that this equation is linear in the components of $H$ ! For each $i$, we have three linear equations for the components of $H$. However, it turns out that one of those three equations is redundant, so each $i$ really only gives us two equations. All together, we have eight linear equations.
$H$ is found by finding a nonzero solution of this system of eight equations. Note that this only determines $H$ up to a scalar multiple. However, that's to be expected, because in homogeneous coordinates, any nonzero scalar multiple of $H$ represents the same projective transformation as $H$.
Best Answer
Here's a way of doing what I think you want to do.
First, find two perpendicular unit-vectors that are parallel to the plane (i.e. perpendicular to the normal vector). You might want to use Gram-Schmidt here. Call the vectors $u_1$ and $u_2$. These are the "axes" along the plane you're projecting onto.
Any point on the plane can be written in the form $a_1u_1 + a_2u_2$ for some coordinates $(a_1,a_2)$. To find the coordinates of the orthogonal projection of a vector $x$ onto your plane, calculate $$ a_1 = u_1^Tx\\ a_2 = u_2^Tx $$ Or, all together: $$ \pmatrix{a_1\\a_2} = \pmatrix{u_1^T\\u_2^T} \pmatrix{x_1\\x_2\\x_3} $$