Say you have two normal vectors, \vec{n}_1 = (a_1, b_1, c_1), \vec{n}_2 = (a_2, b_2, c_2) as well as two position vectors \vec{p}_1 = (x_1, y_1, z_1), \vec{p}_2 = (x_2, y_2, z_2). The intersection of two planes defined by these vectors is all position vectors \vec{x} = (x, y, z) that satisfy
\left\{\begin{array}{ccc}<br />
(\vec{x} - \vec{p}_1) \cdot \vec{n}_1 = 0\\<br />
(\vec{x} - \vec{p}_2) \cdot \vec{n}_2 = 0\end{array}\right..<br />This is equivalent to solving the linear algebra problem\left [\begin{array}{ccc}<br />
a_1 & b_1 & c_1 \\<br />
a_2 & b_2 & c_2 \end{array}\right] \left [\begin{array}{ccc}x \\ y \\ z\end{array}\right] = \left [\begin{array}{ccc}<br />
\vec{p}_1 \cdot \vec{n_1} \\<br />
\vec{p}_2 \cdot \vec{n_2} \end{array}\right] which can be written as \vec A \vec x = \vec b.
Assuming a solution exists (Rouché-Capelli theorem), the solution should be in the form \vec x = k \vec r + \vec p. Hence \vec A (k \vec r + \vec p) = \vec b where \vec p is a particular solution and \vec r is a homogeneous solution that spans the null space of \vec A. As has been pointed out, \vec r = \vec{n}_1 \times \vec{n}_2, \vec{r} \neq \vec 0 is one such solution. In the case when the two normal vectors are parallel, then any non-zero vector spanned by the basis of one of the planes will do like say a vector projection of one of the standard unit vectors onto the plane (keep switching to find one that's non-zero).
To find a particular solution, simply remove a linearly dependent column vector and set it's corresponding variable to zero and solve the matrix equation. For instance, if the first column is linearly dependent, remove it and let x = 0. You then solve the equation
\left [\begin{array}{ccc}<br />
b_1 & c_1 \\<br />
b_2 & c_2 \end{array}\right] \left [\begin{array}{ccc}y \\ z\end{array}\right] = \left [\begin{array}{ccc}<br />
\vec{p}_1 \cdot \vec{n_1} \\<br />
\vec{p}_2 \cdot \vec{n_2} \end{array}\right] to get a particular solution.
EDIT:
Alternatively, you can get everything by reducing the augmented matrix [\vec A|\vec b] to a reduced row echelon form. Then, each leading coefficient is located on a column i. The ith coordinate of the particular solution is therefore the augmented element on the same row. The coordinates of the homogeneous solution is given by the negatively inverted value of the element on the same row in a linearly dependent column j, with the jth coordinate being 1. For instance, the augmented matrix
\left [\begin{array}{ccc}<br />
1 & 2 & 3 & 5\\<br />
1 & 2 & 4 & 9\end{array}\right]has the reduced row echelon form
\left [\begin{array}{ccc}<br />
1 & 2 & 0 & -7\\<br />
0 & 0 & 1 & 4\end{array}\right] hence the full solution is
\vec x = \left [\begin{array}{ccc}<br />
-7\\0\\4\end{array}\right] + k<br />
\left [\begin{array}{ccc}<br />
-2\\1\\0\end{array}\right]<br />
