Pxy[u_, ax_, bx_, cx_, dx_, ay_, by_, cy_, 
        dy_] :=  {ax*u^3 + bx*u^2 + cx*u  + dx,  
        ay*u^3 + by*u^2 + cy*u + dy}

P0 = {x0, y0};
P1  =  {x1, y1};
P2  =  {x2, y2};
P3  =  {x3, y3};
P00 = {x00, y00};
P33  =  {x33, y33};

v  =  
        Solve[
{Pxy[0, ax, bx,  cx, dx, ay, by, cy, 
                dy]  ==  
              P0, 
Pxy[1/3, ax, bx,  cx, dx, ay, by, cy, 
                dy] == P1, 
 
            Pxy[2/3, ax, bx,  cx, dx, ay, by, cy, dy] == 
              P2, 
Pxy[1, ax, bx,  cx, dx, ay, by, cy, 
                dy] == P3}, 
{ax, bx,  cx, dx, ay, 
            by, cy, dy}];

v

w  =  
        Solve[
{Pxy[0, ax, bx,  cx, dx, ay, by, cy, 
                dy]  ==  
              P0, 
(D[
                    Pxy[u, ax, bx,  cx, dx, ay, by, cy, dy], u] /. 
                  u -> 0) == 
              P00, 
(D[ 
                    Pxy[u, ax, bx,  cx, dx, ay, by, cy, dy], u] /. 
                  u -> 1) == 
              P33, 
Pxy[1, ax, bx,  cx, dx, ay, by, cy, 
                dy] == P3}, 
{ax, bx,  cx, dx, ay, 
            by, cy, dy}];

w

Bezier[u_] := 
      For[k = 0; q = 0, k ≤ 3, k++, 
        q =  q + 
            Binomial[3, k] ((1 - u)^(3 - k))*(u^k)*
              Q[[k + 1]]]

Q = {Q1, Q2, Q3, Q4}

Bezier[1/3]

q

Bezier[0.001]

q

Bezier[0.9999]

q