Simple Intersection Algorithm for Ray and 3D Bézier Curves of Varying Thickness
You say that you don't need the fastest possible algorithm, so instead of solving it analytically, you might try considering the numeric alternatives, which in this case could mean ray marching.

How you could do this is take steps down the ray and test if the point is inside or outside of this bezier tube.

You would iterate down the ray until you either find an intersection or leave the bounding box of this shape (or maybe when you take a maximum number of steps, or some other fail case).

The challenge then only becomes "how can I tell if a point is inside or outside of the shape?" which is a lot simpler.  If you had constant thickness you would just find the closest point on the bezier to your point and do a distance comparison.

Some details you might be interested in:


taking smaller steps, and more of them, will give you more accurate results, but will be more expensive computationally. That makes this algorithm tuneable for quality vs speed.
if the step size is too large or the tube is too thin, you may completely miss the tube when there is really a hit.
when/if you do find a hit, you are finding the first point inside the shape. The real answer is somewhere between the last sample and the current sample.  You could refine the hit location with a binary search or some other method if you wanted to.
you can likely increase the step size as you move farther from the camera. Smaller details are harder to see at greater distances.
you might be able to calculate or approximate the gradient of your point vs bezier tube function as you ray march. If so, this gradient can be used to give you a lower bound on the distance from your point to the bezier tube, which would let you take larger steps through empty space.


Regarding the "secondary ray" issue you mention like for shadow rays or reflection rays that comes up with both analytical and numerical solutions.  Getting a (probably) less precise answer than a numerical solution does make the problem worse, but there should be ways to get around it.

For instance, after you detect a hit, you could find a point on the surface near the hit even if it isn't where the hit was (like, push the hit point to the surface, in whatever direction is easiest to compute), add an epsilon (needed for analytical solutions anyways) and spawn secondary rays from there.

Another way could be to have every object have an ID and when ray marching / tracing, don't count intersections with the same ID as what you just hit.  This only really works for opaque convex objects though.