Here are 6 lines of code (excluding the import and using statements) that define linear, quadratic and cubic Bezier splines using higher order functions.
using Plots
import Base.+, Base.*, Base./
+(f, g) = x -> f(x) + g(x)
*(t::Number, g) = x -> t * g(x)
interpolate(a, b) = t -> (1.0-t)*a + t*b
b1(p1, p2) = interpolate(p1, p2)
b2(p1, p2, p3) = interpolate(b1(p1, p2), b1(p2, p3))
b3(p1, p2, p3, p4) = interpolate(b2(p1, p2, p3), b2(p2, p3, p4))
The idea here is that a Bezier spline of order n is defined by n+1 reference points as an interpolation between the two splines of order n-1 that are defined respectively with the first n points or the last n points. The code renders this definition literally and thus is easy to match up to the intuitive definition.
More importantly, the machine-level code produced by using this very abstract definition is just as compact and efficient as the code produced by an explicit polynomial implementation of a cubic spline. In fact, you can even derive that polynomial definition by giving symbolic values to the spline instead of numerical values. If you do that, the result is not a numerical position, but is the polynomial expressing the spline. You can mix this up by providing numerical values for the points defining the spline, but giving a symbolic value for the interpolating variable t.
As a usage note, to compute the position along a cubic spline defined by points p1, p2, p3, and p4, use b3(p1, p2, p3, p4)(t)(t)(t) where t ranges from 0 to 1. The first application gives you an interpolated quadratic spline. When you feed t to that, you get an interpolated line. When you feed t to that line, you get a point which is your answer. You can use this multiple application to produce pretty visualizations that explain how these splines are constructed, but you mostly would probably rather hide the repeated application in a convenience function.