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.
You probably want to adjust your intuition a bit: the intersection can be just a point, and indeed this is the typical case.
In $\mathbb R^4$, a two-dimensional plane represents two linear (technically “affine”) constraints. So intersecting two planes results in four constraints, which (assuming the constraints are independent) yields a single point.
Contrast to the situation in $\mathbb R^3$ where a plane only determines 1 constraint, so two planes determine 2 constraints, leaving enough freedom (in the generic case) to define a line.
So it’s your 3D intuition tricking you into thinking that two planes intersect in a line. It’s similar to how two lines in $\mathbb R^2$ intersect in a point, but in $\mathbb R^3$ the usual case is that two lines don’t intersect at all.
Best Answer
Further to my comment above.
To find the intersection of two planes defined by the equations $ax+by+cz=d$ and $a'x+b'y+c'z=d'$, eliminate two of the variables, using the third, say $z$, as a parameter. You need to solve a $2\times2$ linear system in $x,y$:
$$\left\{\begin{matrix}ax+by&=&d-cz\\ a'x+b'y&=&d'-c'z\end{matrix}\right.$$
This will yield a parametric equation for the line of interection
$$\left\{\begin{matrix}x=\alpha z+\beta\\y=\alpha' z+\beta'\\z=z\end{matrix}\right.$$
It's not always possible to choose $z$ as the parameter (for instance if the two planes have equations $x=0$ and $z=0$, the parameter will be $y$). So in general, write the parametric equation of the line as
$$\left\{\begin{matrix}x&=&\alpha t+\beta\\y&=&\alpha' t+\beta'\\z&=&\alpha'' t+\beta''\end{matrix}\right.$$
Then, to check that this line lies on any other plane $a''x+b''y+c''z=d''$, simply replace $x,y,z$ by $\alpha t+\beta$, $\alpha' t+\beta'$ and $\alpha'' t+\beta''$ respectively, and verify that the equation holds.