diff --git a/src/bSplines/basis.jl b/src/bSplines/basis.jl
index 8e1e67f..3f72af3 100644
--- a/src/bSplines/basis.jl
+++ b/src/bSplines/basis.jl
@@ -120,11 +120,12 @@ function (basis::Union{Bspline,CurrySchoenberg})(evalpoints)
 end
 
 
-struct pAlloc{T<:Real}
+struct pAlloc{T<:Real,F<:Int}
     B::Matrix{T}
     N::Vector{T}
     left::Vector{T}
     right::Vector{T}
+    spanVec::Vector{F}
 end
 
 
@@ -136,7 +137,7 @@ Evaluate B-spline basis at all evalpoints.
 function (basis::Union{Bspline,CurrySchoenberg})(evalpoints, prealloc::pAlloc)
 
     numBasis = numBasisFunctions(basis)
-    knotSpan = findSpan(numBasis, evalpoints, basis.knotVec, basis.degree)
+    knotSpan = findSpan!(prealloc.spanVec, numBasis, evalpoints, basis.knotVec, basis.degree)
 
     return basisFun!(prealloc, knotSpan, evalpoints, basis)
 end
@@ -155,12 +156,27 @@ function preAlloc(degree::Int, uVector::Vector{T}) where {T}
     left  = zeros(T, degree + 1)
     right = zeros(T, degree + 1)
 
-    return pAlloc(B, N, left, right)
+    spanVec = similar(uVector, eltype(degree))
+
+    return pAlloc(B, N, left, right, spanVec)
+end
+
+
+"""
+    findSpan(b::T, uVec, kVec, p::T)
+
+Allocate memory and call findSpan!    
+"""
+function findSpan(b::T, uVec, kVec, p::T) where {T<:Int}
+
+    spanVec = similar(uVec, eltype(b))
+
+    return findSpan!(spanVec, b, uVec, kVec, p)
 end
 
 
 """
-    findSpan(n::Int, u, knotVector)
+    findSpan!(spanVec, b::T, uVec, kVec, p::T)
 
 Find the spans of a B-spline knot vector at the parametric points 'u', where 'b' is the number of basis functions (control points).
 
@@ -170,23 +186,34 @@ Modification of Algorithm A2.1 from 'The NURBS Book' p. 68.
 
 Assumption that the knotVector is open! (the first and last knot are repeated degree + 1 times)
 """
-function findSpan(b::Int, u, knotVector, p)
+function findSpan!(spanVec, b::T, uVec, kVec, p::T) where {T<:Int}
 
-    # check input
-    #(minimum(u) < knotVector[1]) && (maximum(u) > knotVector[end]) && error("Some value is outside the knot span")
+    for (j, u) in enumerate(uVec)
 
-    # initialize the vector containing the indices 
-    spanVec = similar(u, eltype(b))
-
-    for (j, uEval) in enumerate(u)
-
-        # special case: uEval == to last knot vector entry
-        if uEval == knotVector[end]
+        # --- Special case 
+        if u == kVec[end]
             spanVec[j] = b
             continue
         end
 
-        spanVec[j] = findlast(knotVector .≤ uEval)
+        # --- Do binary search 
+        low  = p + 1
+        high = b + 1
+
+        mid = Base.midpoint(low, high)
+
+        while u < kVec[mid] || u >= kVec[mid + 1] # is u in interval [ kVec[mid]...kVec[mid+1] ) ?
+
+            if u < kVec[mid]
+                high = mid
+            else
+                low = mid
+            end
+
+            mid = Base.midpoint(low, high)
+        end
+
+        spanVec[j] = mid
     end
 
     return spanVec
diff --git a/src/bSplines/basisDerivatives.jl b/src/bSplines/basisDerivatives.jl
index d021a97..97dae91 100644
--- a/src/bSplines/basisDerivatives.jl
+++ b/src/bSplines/basisDerivatives.jl
@@ -81,13 +81,14 @@ function (basis::Bspline)(evalpoints, k::Int)
 end
 
 
-struct pAllocDer{T<:Real,L}
+struct pAllocDer{T<:Real,F<:Int,L}
     dersv::Array{T,L}
     ders::Matrix{T}
     ndu::Matrix{T}
     left::Vector{T}
     right::Vector{T}
     a::Matrix{T}
+    spanVec::Vector{F}
 end
 
 
@@ -101,7 +102,7 @@ function (basis::Bspline)(evalpoints::Vector{T}, k::Int, prealloc::pAllocDer) wh
     #basis.divMax < 0 && error("The k-th derivative has to be k ≥ 0!")
 
     numBasis = numBasisFunctions(basis)
-    knotSpan = findSpan(numBasis, evalpoints, basis.knotVec, basis.degree)
+    knotSpan = findSpan!(prealloc.spanVec, numBasis, evalpoints, basis.knotVec, basis.degree)
 
     return derBasisFun!(prealloc, knotSpan, basis.degree, evalpoints, basis.knotVec, k)
 end
@@ -122,7 +123,9 @@ function preAllocDer(degree::Int, evalpoints::Vector{T}, numberDerivatives::Int)
     right = zeros(T, degree + 1)
     a     = zeros(T, 2, degree + 1)
 
-    return pAllocDer(dersv, ders, ndu, left, right, a)
+    spanVec = similar(evalpoints, eltype(degree))
+
+    return pAllocDer(dersv, ders, ndu, left, right, a, spanVec)
 end
 
 
diff --git a/src/nurbs/surfaces.jl b/src/nurbs/surfaces.jl
index b9fce6e..4bc2fe5 100644
--- a/src/nurbs/surfaces.jl
+++ b/src/nurbs/surfaces.jl
@@ -94,9 +94,9 @@ Convenience function to compute points on all k derivatives of a NURBSsurface.
 )
 
 
-struct pAllocNURBSsuface{T<:Real,L}
-    preallocU::pAllocDer{T,L}
-    preallocV::pAllocDer{T,L}
+struct pAllocNURBSsuface{T<:Real,F<:Int,L}
+    preallocU::pAllocDer{T,F,L}
+    preallocV::pAllocDer{T,F,L}
     surfaces::Matrix{Matrix{SVector{3,T}}}
     w::Matrix{Matrix{T}}
 end
@@ -187,12 +187,12 @@ function surfaceDerivativesPoints!(
 
     # u-direction: determine the basis functions evaluated at uVector 
     nbasisFun = length(uKnotVector) - uDegree - 1
-    uSpan = findSpan(nbasisFun, uVector, uKnotVector, uDegree)
+    uSpan = findSpan!(preallocU.spanVec, nbasisFun, uVector, uKnotVector, uDegree)
     Nu = derBasisFun!(preallocU, uSpan, uDegree, uVector, uKnotVector, k)
 
     # v-direction: determine the basis functions evaluated at vVector
     nbasisFun = length(vKnotVector) - vDegree - 1
-    vSpan = findSpan(nbasisFun, vVector, vKnotVector, vDegree)
+    vSpan = findSpan!(preallocV.spanVec, nbasisFun, vVector, vKnotVector, vDegree)
     Nv = derBasisFun!(preallocV, vSpan, vDegree, vVector, vKnotVector, k)
 
     # initialize