Obtaining the normals for a point on the Bezier curve is actually quite straightforward, since the normals are simply perpendicular to the tangent function (oriented in the plane of the direction of movement for the curve) and the tangent Bezier function, the curve is actually just another Bezier curve, 1 order lower. Let there be a normal for the cubic Bezier curve. The regular function c (a, b, c, d) is the coordinates of the curve in one dimension:
function computeBezier (t,a,b,c,d) { return a * (1-t)³ + 3 * b * (1-t)² * t + 3 * c * (1-t) * t² + d * t³ }
Note that the Bezier curves are symmetrical, the only difference between t vs 1-t is that the end of the curve represents the “beginning”. Using a * (1-t)³ means that the curve starts with a . Using a * ³t , you start with d .
So, let's define a fast curve with the following coordinates:
a = (-100,100,0) b = (200,-100,100) c = (0,100,-500) d = (-100,-100,100)

To get the normal value for this function, we first need derivative :
function computeBezierDerivative (t,a,b,c,d) { a = 3*(b−a) b = 3*(cb) c = 3*(dc) return a * (1-t)² + 2 * b * (1-t) * t + 3 * c * t² }
Done. The derivative calculation is stupidly simple (a fantastic property of Bezier curves).
Now, to get the normal, we need to take the normalized tangent vector for some value of t and rotate it by a quarter of a revolution. We can turn it into several directions, so the further restriction is that we want to turn it only into a plane, which is determined by the tangent vector, and the tangent vector "next to it", an infinitesimal interval.
The tangent vector for any Bezier curve is formed simply due to the fact that you have multidimensional measurements and their evaluation separately, so for a 3D curve:
| computeBezierDerivative(t, x values) | |x'| Tangent(t) = | computeBezierDerivative(t, y values) | => |y'| | computeBezierDerivative(t, z values) | |z'|
Again, simple enough to calculate. To normalize this vector (or virtually any vector), we simply perform vector division along the length:
|x'| NormalTangent(t) = |y'| divided by sqrt(x'² + y'² + z'²) |z'|
So let's paint in green:

The only trick is now to find a plane in which you can rotate the tangent vector to turn the tangent into normal. We know that we can use another value of t, arbitrarily close to the one we want, and turn it into a second tangent vector located near one point to find a plane with arbitrary correctness, so we do the following:
Given the starting point f(t1)=p , we take the point f(t2)=q with t2=t1+e , where e is some small value, such as 0.001 - this point q has the derivative q' = pointDerivative(t2) , and in order to do things, it’s easier for us to transfer this tangent vector to the smaller bit on pq so that both vectors “start” with p . Pretty simple.
Now we have two vectors extending in the same coordinate: our real tangent and the "closest" tangent of a point that is so close that it can be the same point. Fortunately, because of the way the Bezier curves work, this second tangent is never the same, but slightly different, and “slightly different” is all we need: if we have two normalized vectors, starting with one at the same point, but pointing in different directions, we can find the axis along which we need to rotate it to get the other, simply using the cross product between them, and therefore we can find that they both pass.
Ordering questions: we calculate tangent₂ × tangent₁ because if we calculate the tangent ₁ × tangent₂ we will calculate the rotation axis and the resulting normals in the "wrong" direction, Correction, it is literally just "take the vector, multiply by -1" at the end but why fix it after we can fix it, here. Let's look at these rotation axes in blue:

Now we have everything we need: to turn our normalized tangent vectors into normal vectors, you just need to rotate them around the axes that we just found a quarter turn. If we turn them in one direction, we get the normals, if we change them to the other, we get the normals of the return trip.
For arbitrary rotation around the axis in 3D, this work is probably laborious, but not complicated , and quarter turns, as a rule, that they greatly simplify the mathematics: to rotate a point over our axis of rotation c , the rotation matrix is equal to:
| c₁² c₁*c₂ - c₃ c₁*c₃ + c₂ | R = | c₁*c₂ + c₃ c₂² c₂*c₃ - c₁ | | c₁*c₃ - c₂ c₂*c₃ + c₁ c₃² |
Where the indices 1, 2, and 3 are really only the x, y, and z components of our vector. So it’s still easy, and all that’s left is the rotation matrix of our normalized tangent:
n = R * Tangent "T"
What is:
| T₁ * R₁₁ + T₂ * R₁₂ + T₃ * R₁₃ | |nx| n = | T₁ * R₂₁ + T₂ * R₂₂ + T₃ * R₂₃ | => |ny| | T₁ * R₃₁ + T₂ * R₃₂ + T₃ * R₃₃ | |nz|
And we have the normal vector (s) we need. Excellent!
Except for what we can do better : in the same way that the vector c was perpendicular to two tangents, our normal n is perpendicular to both c and regular tangency, so we can use the transverse product a second time to find the normal:
|nx| n = c × tangent₁ => |ny| |nz|
This will give us exactly the same vector, with less work.

And if we want internal normals, this is the same vector, just multiply by -1:

Pretty easy as soon as you learn the tricks! And finally, because the code is always useful, this meaning is Processing I used to make sure that I am telling the truth.