(2D) How to make Bezier curve line has width by using python?

If I have a Bezier curve, how could I make it have width, and how to get vertices of its contour?

My attempt: I have plot a Bezier curve by using one start point, two control points and one end points, their coordinates are:

p0_ = [11, -0.45]
p1_ = [13.5, -0.45]
p2_ = [13.5, -4]
p3_ = [16, -4]

the figure is as shown below:

But I want to make it have width, lets say width=0.7, then I change the original coordinates into:

width = 0.7
p01_ = [11, -0.45 + width / 2]
p11_ = [13.5, -0.45 + width / 2]
p21_ = [13.5, -4 + width / 2]
p31_ = [16, -4 + width / 2]
p02_ = [11, -0.45 - width / 2]
p12_ = [13.5, -0.45 - width / 2]
p22_ = [13.5, -4 - width / 2]
p32_ = [16, -4 - width / 2]

And then I plot new figure, I find it looks quite strange. It is not uniform.

Obviously, this is not the one I want, I guess the coordinates are wrong, But I do not know how to get correct one? What I want is like the one shown below:

You see, everywhere is uniform.

This follows the prescription in MBo's comment.

For each point on the Bezier curve you compute the normal and then go width/2 along that normal to each side of the Bezier curve.

To compute the normals you could take successive line segments along your discretised Bezier curve. However, I think you will get a smoother answer (especially at the end points of each Bezier) if you obtain the normals by differentiating the original Bezier curve r(t) to get a tangent dr/dt and then get a normal as k cross dr/dt, where k is a unit vector in the z direction. This is what is done below.

If you want a filled or hatched region then take the left (xL,yL) and right (xR,yR) vertices and create a polygon, which you can then fill. (You will have to reverse either the left or right vertices so that you continue to traverse the polygon in the same sense.)

import numpy as np
import matplotlib.pyplot as plt

#-----------------------------

class Bezier:
    '''class for Cubic Bezier curve'''
    def __init__( self, q0, q1, q2, q3 ):
        # control points
        self.p0 = np.array( q0 )
        self.p1 = np.array( q1 )
        self.p2 = np.array( q2 )
        self.p3 = np.array( q3 )
        # normals (z cross p)
        self.n0 = np.array( [-q0[1], q0[0] ] )
        self.n1 = np.array( [-q1[1], q1[0] ] )
        self.n2 = np.array( [-q2[1], q2[0] ] )
        self.n3 = np.array( [-q3[1], q3[0] ] )

    def pt( self, t ):
        return (1-t)**3 * self.p0 + 3*t*(1-t)**2 * self.p1 + 3*t**2*(1-t) * self.p2 + t**3 * self.p3

    def sides( self, t, width ):
        normal = -3*(1-t)**2 * self.n0 + (3*(1-t)**2-6*t*(1-t)) * self.n1 + (6*t*(1-t)-3*t**2) * self.n2 + 3*t**2 * self.n3
        normal /= np.linalg.norm( normal )
        point = self.pt( t )
        return point + ( width / 2 ) * normal, point - ( width / 2 ) * normal

#-----------------------------

def getCurves( BezierData, n, width ):
    x = [];   y = [];   xL = [];   yL = [];   xR = [];   yR = []
    for B in BezierData:
        arc = Bezier( B[0], B[1], B[2], B[3] )
        for t in np.linspace( 0.0, 1.0, n ):
            P = arc.pt( t );   x.append( P[0] );   y.append( P[1] )
            L, R = arc.sides( t, width );   xL.append( L[0] );   yL.append( L[1] );   xR.append( R[0] );   yR.append( R[1] )
    return x, y, xL, yL, xR, yR
    
#-----------------------------

data = [
        [ [11,-0.45], [13.5, -0.45], [13.5, -4], [16, -4] ]
       ]
n = 100           # points per Bezier
width = 0.7       # width
x, y, xL, yL, xR, yR =  getCurves( data, n, width )
plt.plot( x , y , 'k:' )
plt.plot( xL, yL, 'b-' )
plt.plot( xR, yR, 'b-' )
plt.show()