There is some fairly obscure theory that you could use to shorten your work. Given a polynomial $f$ of degree $n$, there is a unique symmetric multi-affine function $F$ of $n$ variables such that $F(t,t,\ldots, t) = f(t)$ for all $t$. The function $F$ is called the polar form or the blossom of $f$. This concept is well-known among people who work with Bezier curves, and also in classical algebraic geometry. More info here.
Some examples of blossoms are:
$$f(t) = t \Longrightarrow F(u,v,w) = \tfrac13(u+v+w)$$
$$f(t) = t^2 \Longrightarrow F(u,v,w) = \tfrac13(vw + wu + uv)$$
$$f(t) = (1-t)^3 \Longrightarrow F(u,v,w) = (1-u)(1-v)(1-w)$$
In particular, the polar form of your cubic curve is
$$
F(u,v,w) = \Bigl(\; \tfrac13(u+v+w),\;
\tfrac13(vw + wu + uv),\; uvw \;\Bigr)
$$
Polar forms have some nice geometric properties.
For example, for a parametric quadratic curve, the polar form $F(u,v)$ gives the intersection point of the tangents at the parameter values $u$ and $v$.
Similarly, for a non-planar parametric cubic curve, the polar form $F(u,v, w)$ gives the intersection point of the osculating planes at the parameter values $u$, $v$, and $w$. I see that you already rediscovered this for yourself (after a lot of algebra, I would guess). To be honest, I can't remember where I learned this, and I can't find a proof.
To show that the four points are coplanar, you can look at the determinant
$$\left|\matrix{
\tfrac13(u+v+w) & \tfrac13(vw + wu + uv) & uvw & 1 \\
u & u^2 & u^3 & 1 \\
v & v^2 & v^3 & 1 \\
w & w^2 & w^3 & 1
}\right|
$$
This determinant gives the volume of the tetrahedron defined by the four points, so it is zero precisely when the four points are coplanar. There is probably some clever way to show that the determinant is zero, but I don't see it, right now. I expect that it's easy to show by brute-force algebra or by using a system like Mathematica.
Alternatively (and much better, I think) …
Let $P(u)$, $P(v)$, and $P(w)$ be the three points on the curve. Then a formula on page 16 of Ramshaw’s document says that
$$
F(u,v,w) = \frac{(w-v)^2}{3(w-u)(u-v)}P(u) +
\frac{(w-u)^2}{3(w-v)(v-u)}P(v) +
\frac{(v-u)^2}{3(v-w)(w-u)}P(w)
$$
It’s straightforward to show that the three coefficients in this formula add up to 1. So, the formula says that the intersection point $F(u,v,w)$ is an affine combination of the the three curve points, and therefore lies in the plane containing them.
Note that this last formula is valid for any parametric cubic curve, not just the very simple one that you’re considering. Pretty neat, huh. We can thank Monsieur de Casteljau for this.
Best Answer
Just apply Serret-Frenet.
WLOG the common point is the origin and the curve is parametrised by arclength $s$. The given condition says $\mathbf{r}\cdot(\mathbf{r}'\times\mathbf{r}'')=0$.
Since $\mathbf{r}'\times\mathbf{r}''=\mathbf{T}\times\mathbf{T}'=\kappa\mathbf{T}\times\mathbf{N}=\kappa\mathbf{B}$, either we have $\kappa=0$ on an interval (so we are done), or $\mathbf{r}\cdot\mathbf{B}=0$ on an interval.
Differentiating $\mathbf{r}\cdot\mathbf{B}=0$ gives $\mathbf{r}\cdot\tau\mathbf{N}=0$ so again either zero torsion on an interval (so we are done) or $\mathbf{r}\cdot\mathbf{N}=0$ on an interval. Since $\mathbf{r}\cdot\mathbf{N}=\mathbf{r}\cdot\mathbf{B}=0$ we have $\mathbf{r}=(\mathbf{r}\cdot\mathbf{T})\mathbf{T}$ so $\kappa=0$.