The general equation of a plane in 3-D is given by $$\mathbf{(p-p_0).n}=0$$
where $\mathbf{p}$ is any general point on the plane, and $\mathbf{p_0}$ is any known point on the plane. $\mathbf{n}$ is a vector normal to the plane.
The equation of a line is given by $$\mathbf{p = l_0+}t\mathbf{l} $$
where $\mathbf{l_0}$ is any point on the line.
If the line lies in the plane, it must satisfy two conditions-
It must be perpendicular to the normal to the plane i.e. $\mathbf{l.n}=0$
$\mathbf{l_0}$ must lie in the plane i.e. satisfy the plane's equation. So, $\mathbf{(l_0-p_0).n} = 0$
You can calculate $\mathbf{p_0,l_0}$ very easily with the information you have.
$\mathbf{p_0}$ can be calculated by choosing any two of $x,y,z$ and finding the third to satisfy the equation of the plane. $\mathbf{l_0}$ is precisely $\vec O$ that you already know. You can use them to cross-check whether the fit is good or not. But they do not depend on $\mathbf{l}$. So, they are secondary, and may be used as a sanity check later on.
Suppose you are given the equation of the line as $ax+by+cx+d = 0$. You can recast it as $[(x,y,z)-\mathbf{p_0}].(a,b,c)$. So, $(a,b,c)$ is your normal vector.
Suppose you have $m$ planes with normals $\mathbf{n_1,n_2,\ldots,n_m}$. Then, your overall constraints are $$\mathbf{l.n_1} = 0 \\ \mathbf{l.n_2} = 0 \\ \ldots \\\mathbf{l.n_m}=0$$.
It is a system of linear equations with 3 variables and $m$ equations. You need to find out $\mathbf{l}$, which can be found up to a constant factor by the standard least squares method, provided $m>3$, which should not be a worry for you. If you code carefully enough, you can implement all of it in matrix terms.
Hope it helps.
Choose orthonormal bases $(E_1, E_2)$ and $(F_1, F_2)$ for the two planes. Then, the (hyper)volume of the parallelepiped spanned by the $E_1, E_2, F_1, F_2$ on the one hand is $\sin \theta$ where $\theta$ is the angle between the two planes, and on the other hand is (the absolute value of) the determinant of the matrix given by adjoining the four vectors:
$|\det[E_1 \, E_2 \, F_1 \, F_2]|$.
This quantity is independent of the choices of orthogonal bases, and in fact, we can even take the bases of the planes to be any that span parallelograms of unit area, so we don't need to produce an orthogonal basis.
If we start out with any bases $(E_i)$ or $(F_i)$, say ones that span areas $\lambda$ and $\mu$, we can always normalize them by rescaling vectors, but we might as well build this into our equation: The bases $(\lambda^{-1} E_1, E_2)$ and $(\mu^{-1} F_1, F_2)$ both span parallelograms of unit area, so for general bases the hypervolume is
$\sin \theta = |\det[\lambda^{-1} E_1 \, E_2 \, \mu^{-1} F_1 \, F_2]| = \frac{|\det[E_1 \, E_2 \, F_1 \, F_2]|}{\lambda \mu}$,
which we might write as
$\color{red}{\sin \theta = \dfrac{|\det[E_1 \, E_2 \, F_1 \, F_2]|}{|E_1 \wedge E_2| |F_1 \wedge F_2|}}$,
where $|G_1 \wedge G_2|$ denotes the area of the parallelograms spanned by $G_1, G_2$.
This generalizes readily to formulas the angle between $k$-planes and $(n - k)$-planes in vector spaces of dimension $n$ (try this for the familiar situation $n = 2$, $k = 1$), and with just a little more work to finding the angle between $k$- and $l$- planes in vector spaces of dimension more than $k + l$.
Remark This leaves the matter of computing explicitly the areas $|G_1 \wedge G_2|$ of the parallelograms the bases defime. The area $A$ of the parallelogram defined by vectors $H_1 = (x_1, y_1), H_2 = (x_2, y_2)$ in the plane is
$A = \left\vert\det \left(\begin{array}{cc} x_1 & x_2 \\ y_1 & y_2\end{array}\right)\right\vert = |x_1 y_2 - x_2 y_1|$, and its square is $A^2 = (x_1 y_2 - x_2 y_1)^2$, which we can rewrite as
$A^2 = [(x, y) \cdot (x, y)] [(x', y') \cdot (x', y')] - [(x, y) \cdot (x', y')]^2 = (H_1 \cdot H_1) (H_2 \cdot H_2) - (H_1 \cdot H_2)^2$.
Now, the formula $A^2 = (H_1 \cdot H_1) (H_2 \cdot H_2) - (H_1 \cdot H_2)^2$ doesn't depend on coordinates, it just uses the Euclidean structure (namely the dot product $\cdot\,$), so it works just as well for computing the areas of the parallelograms in our original problem, that is, we may write
$|G_1 \wedge G_2|^2 = (G_1 \cdot G_1) (G_2 \cdot G_2) - (G_1 \cdot G_2)^2$
and then take square roots if we like.
Best Answer
Suppose that $\left|\begin{matrix} A_1&B_1 \\ A_2 & B_2\end{matrix}\right| = A_1B_2-B_1A_2\neq 0$. Then you may reformulate the problem as follows:
$$\begin{pmatrix} A_1&B_1 \\ A_2 & B_2\end{pmatrix} \begin{pmatrix} x \\ y\end{pmatrix} = -\begin{pmatrix} C_1 z + D_1 \\ C_2 z+D_2\end{pmatrix} $$ and solve for $x$ and $y$: $$ \begin{pmatrix} x \\ y\end{pmatrix} = -\begin{pmatrix} C'_1 z + D'_1 \\ C'_2 z+D'_2\end{pmatrix} $$ This shows that for any $z=t\in{\Bbb R}$ you get a unique solution for $x$ and $y$. What happens here is that the intersection of the two planes $P_1,P_2$ with the plane $z-t=0$ provides two non-parallel lines (because of the non-zero A-B determinant) in the $x-y$ plane. These two lines therefore have a unique intersection point.
Now, when your A-B determinant above is zero (so your two lines in the $x-y$ plane are parallel) then you may look for a non-zero $B-C$ matrix (and solve for $y,z$) or a non-zero $C-A$ matrix (and solve for $z,x$). If all these determinants are zero then your two original planes are in fact parallel so either the intersection is empty or it is a plane.
Note that the three determinants you compute are in fact the component of the cross-product of normal vectors for the planes, so the cross-product being non-vanishing is indeed a condition for the intersection to be a line.