You just need to fully compile your function: fullycompiledBSplineSurf = Hold@Compile[{{ctrlnets, _Real, 3}, {deg1, _Integer}, {deg2, _Integer}, {knots1, _Real, 1}, {knots2, _Real, 1}, {u, _Real}, {v, _Real}}, Module[{i, j, validnets, row, col}, i = searchSpan[{deg1, knots1}, u]; j = searchSpan[{deg2, knots2}, v]; validnets = Take[ctrlnets, {i - deg1 + 1, i + 1}, {j - deg2 + 1, j + 1}]; row = optimizedNonzeroBasis[deg1, knots1, u]; col = optimizedNonzeroBasis[deg2, knots2, v]; row.Transpose[validnets, {1, 3, 2}].col], CompilationOptions -> {"InlineExternalDefinitions" -> True}, "RuntimeOptions" -> {"EvaluateSymbolically" -> False}, CompilationTarget -> C, RuntimeOptions -> "Speed"] /. DownValues@searchSpan // ReleaseHold; ParametricPlot3D[ fullycompiledBSplineSurf[pts, 3, 3, k1, k2, u, v], {u, 0, 1}, {v, 0, 1}] // AbsoluteTiming [![enter image description here][1]][1] I think I don't need to make any explanation, because all the techniques above have been used in answering your previous questions. [1]: https://i.sstatic.net/x3WOH.png