Let's look at an approach that works for all cubic curves:
$$\vec{p}(t) = \vec{p}_0 + t \vec{p}_1 + t^2 \vec{p}_2 + t^3 \vec{p}_3 \tag{1}\label{1}$$
with $0 \le t \le 1$, $t \in \mathbb{R}$.
The straight line is from $ \vec{p}(0) = \vec{p}_0 $ to $ \vec{p}(1) = \vec{p}_0 + \vec{p}_1 + \vec{p}_2 + \vec{p}_3 $.
Let the direction vector for the line be $(x_E , y_E , z_E ) = \vec{p}(1) - \vec{p}(0) = \vec{p}_1 + \vec{p}_2 + \vec{p}_3$.
First, we translate and rotate the coordinate system so that this line will be on the $x$ axis. The translation is by $-\vec{p}_0$, and the orthonormal rotation matrix is
$$\mathbf{R} = \left [ \begin{array}{ccc}
r_{11} & r_{12} & r_{13} \\
r_{21} & r_{22} & r_{23} \\
r_{31} & r_{32} & r_{33} \end{array} \right ]$$
where$$r_{11} = \frac{ x_E }{\sqrt{ x_E^2 + y_E^2 + z_E^2 }}$$
$$r_{12} = \frac{ y_E }{\sqrt{ x_E^2 + y_E^2 + z_E^2 }}$$
$$r_{13} = \frac{ z_E }{\sqrt{ x_E^2 + y_E^2 + z_E^2 }}$$
We have an underdetermined situation: we do not care how the other two axes end up, as long as the matrix itself is valid (orthonormal). So, use the Gram-Schmidt process to find the second row vector ($r_{21}$, $r_{22}$, and $r_{23}$) starting with any initial vector not parallel to $(x_E , y_E , z_E )$; remember to normalize the result to unit length. The third row vector is the cross product of the two first row vectors.
Translate and rotate the control points, to get new ones:
$$\begin{cases}
\vec{c}_0 = (0, 0, 0) \\
\vec{c}_1 = \mathbf{R} \left ( \vec{p}_1 - \vec{p}_0 \right ) \\
\vec{c}_2 = \mathbf{R} \left ( \vec{p}_2 - \vec{p}_0 \right ) \\
\vec{c}_3 = \mathbf{R} \left ( \vec{p}_3 - \vec{p}_0 \right ) \end{cases}$$
In the transformed coordinates, the curve is now
$$\vec{c} = t \vec{c}_1 + t^2 \vec{c}_2 + t^3 \vec{c}_3 \tag{2}\label{2}$$
Using Cartesian coordinates, we can write $\eqref{2}$ as
$$\begin{cases}
x(t) = t x_1 + t^2 x_2 + t^3 x_3 \\
y(t) = t y_1 + t^2 y_2 + t^3 y_3 \\
z(t) = t z_1 + t^2 z_2 + t^3 z_3 \end{cases}$$
where $\vec{c}_1 = ( x_1 , y_1 , z_1 )$, $\vec{c}_2 = ( x_2 , y_2 , z_2 )$, and $\vec{c}_3 = ( x_3 , y_3 , z_3 )$.
Remember, $0 \le t \le 1$. We also know that $\vec{c}(0) = (0, 0, 0)$, and that $\vec{c}(1) = (L, 0, 0)$, where $L = \sqrt{x_E^2 + y_E^2 + z_E^2}$. This means that $x(t)$ is the distance from origin along the line.
The distance to from point $\vec{c}(t)$ to the line ($x$-axis) is
$$r(t) = \sqrt{ y(t)^2 + z(t)^2 } \tag{3}\label{3}$$
and we are looking for the extrema of $\eqref{3}$ within $0 \le t \le 1$. Because $r(t)$ is a nonnegative function, it reaches its extrema where its square does. Therefore, we can actually look for the extremas of
$$R(t) = r(t)^2 = y(t)^2 + z(t)^2$$
which occur when $t = 0$, $t = 1$, or the derivative is zero;
$$\begin{array}{rl}
\frac{d \, R(t)}{d \, t} = 0 = & t^5 ( 6 y_3^2 + 6 z_3^2 ) + \\
\; & t^4 ( 10 y_2 y_3 + 10 z_2 z_3 ) + \\
\; & t^3 ( 4 y_2^2 + 4 z_2^2 + 8 y_1 y_3 + 8 z_1 z_3 ) + \\
\; & t^2 ( 6 y_1 y_2 + 6 z_1 z_2 ) + \\
\; & t ( 2 y_1^2 + z_1^2 ) \end{array}\tag{4}\label{4}$$
Note that $t = 0$ is always a root; $t = 1$ is mathematically also a root (because of how the coefficients have been defined; but not in general for the above equation!), but rounding errors may perturb the root at $1$ slightly.
$\eqref{4}$ does have analytical solutions, but you'll probably want to use Sage, Maple, Mathematica, or something similar to massage the solutions into more useful form; the solutions I found using Sage are just too long to paste here.
Let's say you find roots $t_i$, $0 \le t_i \le 1$. Pick $t_{max} = t_k$ where $R(t_k) \ge R(t_i)$ for all $i$; i.e. pick the $t$ that yields the largest $R(t)$. Then:
The point where the curve is furthest from the straight line occurs at $t_{max}$.
In original coordinates, the curve is farthest from the straight line at $\vec{p}(t_{max})$.
The distance is then $r(t_{max}) = \sqrt{R(t_{max})}$.
Distance from $\vec{p}(0)$ to the point on the line that is perpendicular to, and closest to the furthest point, is $x(t_{max})$.
