How to create MATLAB code for making bilinear coon's surface patch from cubic Bezier curves

I have following data from which I need to create coon's patch. Matrix of polygon vertices given below:

P0u = \[0,0,0; 0.2357,0.2357,0.3333; 1.1785,0.2357,0.3333; 1.4142,0,0\];

P0w = \[1.4142,0,0; 1.1785,0.2357,0.3333; 1.1785,1.1785,0.3333; 1.4142,1.4142,0\];

P1u = \[1.4142,1.4142,0; 1.1785,1.1785,0.3333; 0.2357,1.1785,0.3333; 0,1.4142,0\];

P1w = \[0,1.4142,0; 0.2357,1.1785,0.3333; 0.2357,0.2357,0.3333; 0,0,0\];

What can I try next?

I have prepared partial MATLAB code but do not know how to proceed further. See the below code:

clc
clear
close all

%% Define input parameters
A1 = [0,0,0];
A2 = [0.2357,0.2357,0.3333];
A3 = [1.1785,0.2357,0.3333];
A4 = [1.4142,0,0];
B1 = [1.4142,0,0];
B2 = [1.1785,0.2357,0.3333];
B3 = [1.1785,1.1785,0.3333];
B4 = [1.4142,1.4142,0];
C1 = [1.4142,1.4142,0];
C2 = [1.1785,1.1785,0.3333];
C3 = [0.2357,1.1785,0.3333];
C4 = [0,1.4142,0];
D1 = [0,1.4142,0];
D2 = [0.2357,1.1785,0.3333];
D3 = [0.2357,0.2357,0.3333];
D4 = [0,0,0];

P0u = [A1;A2;A3;A4]; % Creating matrix for polygon vertices
P0w = [B1;B2;B3;B4]; % Creating matrix for polygon vertices
P1u = [C1;C2;C3;C4]; % Creating matrix for polygon vertices
P1w = [D1;D2;D3;D4]; % Creating matrix for polygon vertices


[r1,s1] = size(P0u); % getting size of matrix in terms of rows and columns
[r2,s2] = size(P0w); % getting size of matrix in terms of rows and columns
[r3,s3] = size(P1u); % getting size of matrix in terms of rows and columns
[r4,s4] = size(P1w); % getting size of matrix in terms of rows and columns

n1 = r1-1; % n+1 represents number of vertices of the polygon
n2 = r2-1; % n+1 represents number of vertices of the polygon
n3 = r3-1; % n+1 represents number of vertices of the polygon
n4 = r4-1; % n+1 represents number of vertices of the polygon

np = 20; % represents number of equi-distance points on the bezier curve
t = linspace(0,1,np);

%% Plot polygon

for k1 = 1:n1
    figure(1)
    plot([P0u(k1,1),P0u(k1+1,1)], [P0u(k1,2),P0u(k1+1,2)], 'r', 'LineWidth', 2)
hold all
grid on
view(45,45)
end

for k2 = 1:n2
    plot([P0w(k2,1),P0w(k2+1,1)], [P0w(k2,2),P0w(k2+1,2)], 'r', 'LineWidth', 2)
end

for k3 = 1:n3
    plot([P1u(k3,1),P1u(k3+1,1)], [P1u(k3,2),P1u(k3+1,2)], 'r', 'LineWidth', 2)
end

for k4 = 1:n4
    plot([P1w(k4,1),P1w(k4+1,1)], [P1w(k4,2),P1w(k4+1,2)], 'r', 'LineWidth', 2)
end

%% Generate the points on the bezier curve

for j1 = 1:np
    P1 = [0,0,0];
    for i1 = 0:n1
        M1(i1+1) = (factorial(n1)/(factorial(i1)*factorial(n1-i1)))*((t(j1))^i1)*((1-t(j1))^(n1-i1));
        P1 = P1 + P0u(i1+1,:)*M1(i1+1);
    end
    Q1(j1,:) = P1;
end

for j2 = 1:np
    P2 = [0,0,0];
    for i2 = 0:n2
        M2(i2+1) = (factorial(n2)/(factorial(i2)*factorial(n2-i2)))*((t(j2))^i2)*((1-t(j2))^(n2-i2));
        P2 = P2 + P0w(i2+1,:)*M2(i2+1);
    end
    Q2(j2,:) = P2;
end

for j3 = 1:np
    P3 = [0,0,0];
    for i3 = 0:n3
        M3(i3+1) = (factorial(n3)/(factorial(i3)*factorial(n3-i3)))*((t(j3))^i3)*((1-t(j3))^(n3-i3));
        P3 = P3 + P1u(i3+1,:)*M3(i3+1);
    end
    Q3(j3,:) = P3;
end

for j4 = 1:np
    P4 = [0,0,0];
    for i4 = 0:n4
        M4(i4+1) = (factorial(n4)/(factorial(i4)*factorial(n4-i4)))*((t(j4))^i4)*((1-t(j4))^(n4-i4));
        P4 = P4 + P1w(i4+1,:)*M4(i4+1);
    end
    Q4(j4,:) = P4;
end

%% Plot the bezier curve from the obtained points

for l1 = 1:np-1
    plot([Q1(l1,1),Q1(l1+1,1)],[Q1(l1,2),Q1(l1+1,2)], 'b', 'LineWidth', 2);
end

for l2 = 1:np-1
    plot([Q2(l2,1),Q2(l2+1,1)],[Q2(l2,2),Q2(l2+1,2)], 'b', 'LineWidth', 2);
end

for l3 = 1:np-1
    plot([Q3(l3,1),Q3(l3+1,1)],[Q3(l3,2),Q3(l3+1,2)], 'b', 'LineWidth', 2);
end

for l4 = 1:np-1
    plot([Q4(l4,1),Q4(l4+1,1)],[Q4(l4,2),Q4(l4+1,2)], 'b', 'LineWidth', 2);
end