Vstep = 1 / (NumberOfVVertices-1);
VstepMatrix = (Vstep³,       0,      0, 0,
                                         0, Vstep²,       0, 0,
                    0,       0, Vstep,  0,
                    0,       0,       0,1);

Ustep = 1 / (NumberOfUVertices-1);
UstepMatrix = (Ustep³,      0,     0, 0,
                    0, Ustep²,     0, 0,
                    0,      0, Ustep, 0,
                    0,      0,     0, 1);

M = (VstepMatrix (B * P * B^T)) * UstepMatrix;

Vtx = M₃₃;
dvVtx = M₃₀ + M₃₁ + M₃₂;
ddvVtx = 6 * M_{30 +} 2 * M₃₁;
dddvVtx = 6 * M₃₀;

duVtx = M₀₃ + M₁₃ + M₂₃;
dvduVtx = M₀₀ + M₀₁ + M₀₂ + M₁₀ + M₁₁ + M₁₂ + M₂₀ + M₂₁ + M₂₂ ;
ddvduVtx = 6 * (M₀₀ + M₁₀ + M₂₀) 2 * (M₀₁ + M₁₁ + M₂₁) ;
dddvduVtx = 6 * (M₀₀ + M₁₀ + M₂₀) ;

dduVtx = 6 * M₀₃ ;
dvdduVtx = 6 * (M₀₀ + M₀₁ + M₀₂) ;
ddvdduVtx = 36 * M₀₀ + 12 * M₀₁ ;
dddvdduVtx = 36 * M₀₀ ;

ddduVtx = 6 * M₀₃ + 2 * M₁₃ ;
dvddduVtx = 6 * (M₀₀ + M₀₁ + M₀₂) + 2 * (M₁₀ + M₁₁ + M₁₂) ;
ddvddduVtx = 36 * M₀₀ + 12 * (M₀₁ + M₁₀) + 4 * M₁₁;
dddvddduVtx = 36 * M₀₀ + 12 * M₁₀;

// TangV

// We don’t forget the factors of vTVector
f3210 = (3, 0, 0, 0,
         0, 2, 0, 0,
         0, 0, 1, 0,
         0, 0, 0, 0);
MM = f3210 * M;

TangentV = MM₃₂;
duTangentV = MM₀₂ + MM₁₂ + MM₂₂;
dduTangentV = MM₀₂ * 6.F + MM₁₂ * 2.F;
ddduTangentV = MM₀₂ * 6.F;

dvTangentV = MM₃₀+MM₃₁;
dvduTangentV = MM₀₀+MM₀₁ + MM₁₀+MM₁₁ + MM₂₀+MM₂₁;
dvdduTangentV = = MM₀₀+MM₀₁ * 6.F + = MM₁₀+MM₁₁ * 2.F;
dvddduTangentV = = MM₀₀+MM₀₁ * 6.F;

ddvTangentV = MM₃₀*2.F;
ddvduTangentV = MM₀₀*2.F + MM₁₀*2.F + MM₂₀*2.F;
ddvdduTangentV = MM₀₀*12.F + MM₁₀*4.F;
ddvddduTangentV = MM₀₀*12.F;

// TangU
// We don’t forget the factors of uTVector
MM = M *f3210;

TangentU = MM₂₃;
duTangentU = MM₀₃ + MM₁₃;
dduTangentU = MM₀₃ * 2.F;

dvTangentU = MM₂₀+MM₂₁+MM₂₂;
dvduTangentU = MM₀₀+MM₀₁+MM₀₂ + MM₁₀+MM₁₁+MM₁₂;
dvdduTangentU = = MM₀₀+MM₀₁+MM₀₂ * 2.F;

ddvTangentU = MM₂₀*6.F+MM₂₁*2.F;
ddvduTangentU = MM₀₀*6.F+MM₀₁*2.F + MM₁₀*6.F+MM₁₁*2.F;
ddvdduTangentU = MM₀₀*12.F+MM₀₁*4.F;

dddvTangentU = MM₂₀*6.F;
dddvduTangentU = MM₀₀*6.F + MM₁₀*6.F;
dddvdduTangentU = MM₀₀ * 12.F;

for(unsigned int vInt=0; vInt < NumberOfVVertices; vInt++)
{

lVtx = Vtx;
lduVtx = duVtx;
ldduVtx = dduVtx;

lTangentV = TangentV;
lduTangentV = duTangentV;
ldduTangentV = dduTangentV;

lTangentU = TangentU;
lduTangentU = duTangentU;

for(unsigned int uInt=0; uInt < NumberOfUVertices; uInt++)
{

output(Vtx);
Normal = Normalize(CrossProduct(TangentV, TangentU));
output(Normal);

lVtx += lduVtx;
lduVtx += ldduVtx;
ldduVtx += ddduVtx;

lTangentV += lduTangentV;
lduTangentV += ldduTangentV;
ldduTangentV += ddduTangentV;

lTangentU += lduTangentU;
lduTangentU += ldduTangentU;

}

// Vertex
Vtx += dvVtx;
dvVtx += ddvVtx;
ddvVtx += dddvVtx;

duVtx += dvduVtx;
dvduVtx += ddvduVtx;
ddvduVtx += dddvduVtx;

dduVtx += dvdduVtx;
dvdduVtx += ddvdduVtx;
ddvdduVtx += dddvdduVtx;

ddduVtx += dvddduVtx;
dvddduVtx += ddvddduVtx;
ddvddduVtx += dddvddduVtx;

//TangV
TangentV += dvTangentV;
duTangentV += dvduTangentV;
dduTangentV += dvdduTangentV;
ddduTangentV += dvddduTangentV;

dvTangentV += ddvTangentV;
dvduTangentV += ddvduTangentV;
dvdduTangentV += ddvdduTangentV;
dvddduTangentV += ddvddduTangentV;

// TangU
TangentU += dvTangentU;
duTangentU += dvduTangentU;
dduTangentU += dvdduTangentU;

dvTangentU += ddvTangentU;
dvduTangentU += ddvduTangentU;
dvdduTangentU += ddvdduTangentU;

ddvTangentU += dddvTangentU;
ddvduTangentU += dddvduTangentU;
ddvdduTangentU += dddvdduTangentU;

}

}