diff --git a/.travis.yml b/.travis.yml
index 1879b29..1f8732f 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -3,6 +3,7 @@ os:
     - linux
 julia:
     - 0.6
+    - 0.7
 before_script:
     - julia --color=yes -e 'Pkg.clone("https://github.com/JuliaFEM/PkgTestSuite.jl.git")'
     - julia --color=yes -e 'using PkgTestSuite; init()'
diff --git a/src/FEMBasis.jl b/src/FEMBasis.jl
index fa38bda..08086ee 100644
--- a/src/FEMBasis.jl
+++ b/src/FEMBasis.jl
@@ -2,10 +2,14 @@
 # License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
 
 module FEMBasis
-#using Logging
+
+using LinearAlgebra
 
 abstract type AbstractBasis end
 
+include("subs.jl")
+include("vandermonde.jl")
+
 include("create_basis.jl")
 export get_reference_element_coordinates, eval_basis!, eval_dbasis!
 
diff --git a/src/create_basis.jl b/src/create_basis.jl
index 0ce35f4..5443930 100644
--- a/src/create_basis.jl
+++ b/src/create_basis.jl
@@ -3,142 +3,97 @@
 
 using Calculus
 
+import Base.parse
+
 function get_reference_element_coordinates end
 function eval_basis! end
 function eval_dbasis! end
 
-function get_dim(p::Expr)
-    pargs = p.args[2:end]
-    :w in pargs && return 3
-    :v in pargs && return 2
-    :u in pargs && return 1
-end
-
-"""
-    function calculate_basis_coefficients(polynomial::String, coordinates::Vararg{Tuple})
-
-Calculate "interpolate coefficient matrix" for some polynomial p.
-
-# Examples
-That is, if we have polynomial p = 1 + u + v and coordinates (0,0), (1,0), (0,1),
-we find A*p such that first row is the first coordinate, second row is second
-coordinate and so on:
-```julia
-julia> p = "1 + u + w"
-julia> X = ((0.0,0.0), (1.0,0.0), (0.0,1.0))
-julia> calculate_basis_coefficient(p, X)
-[1.0 0.0 0.0 # <-- p(0.0,0.0) = 1.0     = [1.0 0.0 0.0] * [1.0, u, v]
- 1.0 1.0 0.0 # <-- p(1.0,0.0) = 1.0 + u = [1.0 1.0 0.0] * [1.0, u, v]
- 1.0 0.0 1.0] # <- p(0.0,1.0) = 1.0 + v = [1.0 0.0 1.0] * [1.0, u, v]
-```
-"""
-function calculate_basis_coefficients(p::Expr, X::Tuple)
-    pargs = p.args[2:end]
-    dim = get_dim(p)
-    nbasis = length(X)
-    sandbox = Module(:__SANDBOX__)
-    vars = (:u, :v, :w)
-
-    # calculate coefficient matrix A and inversion
-    A = zeros(nbasis, nbasis)
-    for i=1:nbasis
-        for j=1:nbasis
-            for (k,c) in enumerate(X[j])
-                code = :( $(vars[k]) = $c )
-                eval(sandbox, code)
-            end
-            A[j,i] = eval(sandbox, pargs[i])
+function calculate_interpolation_polynomials(p, V)
+    basis = []
+    first(p.args) == :+ || error("Use only summation between terms of polynomial")
+    args = p.args[2:end]
+    N = size(V, 1)
+    b = zeros(N)
+    for i in 1:N
+        fill!(b, 0.0)
+        b[i] = 1.0
+        solution = V \ b
+        N = Expr(:call, :+)
+        for (ai, bi) in zip(solution, args)
+            isapprox(ai, 0.0) && continue
+            push!(N.args, simplify( :( $ai * $bi ) ))
         end
+        push!(basis, N)
     end
-    return A 
+    return basis
 end
 
-"""
-"""
-function calculate_interpolation_polynomials(p::Expr, A::Matrix)
-    pargs = p.args[2:end]
-    dim = get_dim(p)
-    nbasis = size(A, 1)
-    invA = inv(A)
-    equations = Expr[]
-    for j=1:nbasis
-        equation = Expr(:call, :+)
-        for i=1:nbasis
-            isapprox(invA[i,j], 0.0) && continue
-            if typeof(pargs[i]) <: Int && pargs[i] == 1
-                push!(equation.args, invA[i,j])
-            else
-                if isapprox(invA[i,j], 1.0)
-                    push!(equation.args, :($(pargs[i])))
-                elseif isapprox(invA[i,j], -1.0)
-                    push!(equation.args, :(-$(pargs[i])))
-                else
-                    push!(equation.args, :($(invA[i,j]) * $(pargs[i])))
-                end
-            end
+function calculate_interpolation_polynomial_derivatives(basis, D)
+    vars = [:u, :v, :w]
+    dbasis = []
+    for N in basis
+        partial_derivatives = []
+        for j in 1:D
+            push!(partial_derivatives, simplify(differentiate(N, vars[j])))
         end
-        push!(equations, equation)
+        push!(dbasis, partial_derivatives)
     end
-    return equations
+    return dbasis
 end
 
-"""
-Calculate derivatives of basis functions with respect to parameters u, v, w.
-"""
-function calculate_interpolation_polynomial_derivatives(basis, dim)
-    equations = Vector[]
-    vars = [:u, :v, :w]
-    nbasis = length(basis)
-    for j=1:nbasis
-        deq = []
-        for k=1:dim
-            #println("∂($(basis[j])) / ∂$(vars[k])")
-            der = simplify(differentiate(basis[j], vars[k]))
-            push!(deq, der)
-        end
-        push!(equations, deq)
-    end
-    return equations
+function create_basis(name, description, X::NTuple{N, NTuple{D, T}}, p::Expr) where {N, D, T}
+    @debug "create basis given antsatz polynomial" name description X p
+    V = vandermonde_matrix(p, X)
+    basis = calculate_interpolation_polynomials(p, V)
+    return create_basis(name, description, X, basis)
 end
 
-function create_basis(name::Symbol, description::String, X::NTuple{nbasis, NTuple{dim, Float64}}, basis, dbasis) where {nbasis,dim}
+function create_basis(name, description, X::NTuple{N, NTuple{D, T}}, basis) where {N, D, T}
+    @debug "create basis given basis functions" name description X basis
+    dbasis = calculate_interpolation_polynomial_derivatives(basis, D)
+    return create_basis(name, description, X, basis, dbasis)
+end
+
+function create_basis(name, description, X::NTuple{N, NTuple{D, T}}, basis, dbasis) where {N, D, T}
+
+    @debug "create basis given basis functions and derivatives" name description X basis dbasis
 
     Q = Expr(:block)
-    for i=1:nbasis
+    for i=1:N
         push!(Q.args, :(N[$i] = $(basis[i])))
     end
 
     V = Expr(:block)
-    for i=1:nbasis
-        for k=1:dim
+    for i=1:N
+        for k=1:D
             push!(V.args, :(dN[$k,$i] = $(dbasis[i][k])))
         end
     end
 
-    # FIXME: this can be done better somehow
-    if dim == 1
+    if D == 1
         unpack = :((u,) = xi)
-    elseif dim == 2
+    elseif D == 2
         unpack = :((u, v) = xi)
     else
         unpack = :((u, v, w) = xi)
     end
 
     code = quote
-        mutable struct $name <: FEMBasis.AbstractBasis
+        struct $name <: FEMBasis.AbstractBasis
         end
 
         function Base.size(::Type{$name})
-            return ($dim, $nbasis)
+            return ($D, $N)
         end
 
         function Base.size(::Type{$name}, j::Int)
-            j == 1 && return $dim
-            j == 2 && return $nbasis
+            j == 1 && return $D
+            j == 2 && return $N
         end
 
         function Base.length(::Type{$name})
-            return $nbasis
+            return $N
         end
 
         function FEMBasis.get_reference_element_coordinates(::Type{$name})
@@ -161,24 +116,3 @@ function create_basis(name::Symbol, description::String, X::NTuple{nbasis, NTupl
     return code
 end
 
-function create_basis(name::Symbol, description::String, X::NTuple{nbasis, NTuple{dim, Float64}}, p::Expr) where {nbasis,dim}
-    A = calculate_basis_coefficients(p, X)
-    basis = calculate_interpolation_polynomials(p, A)
-    dbasis = calculate_interpolation_polynomial_derivatives(basis, dim)
-    return create_basis(name, description, X, basis, dbasis)
-end
-
-function create_basis(name::Symbol, description::String, X::NTuple{nbasis, NTuple{dim, Float64}}, p::String) where {nbasis,dim}
-    create_basis(name, description, X, parse(p))
-end
-
-function create_basis(name::Symbol, description::String, X::NTuple{nbasis, NTuple{dim, Float64}}, basis_::NTuple{nbasis, String}) where {nbasis,dim}
-    basis = [parse(b) for b in basis_]
-    dbasis = calculate_interpolation_polynomial_derivatives(basis, dim)
-    return create_basis(name, description, X, basis, dbasis)
-end
-
-function create_basis(name::Symbol, description::String, X::NTuple{nbasis, NTuple{dim, Float64}}, basis::NTuple{nbasis, Expr}) where {nbasis,dim}
-    dbasis = calculate_interpolation_polynomial_derivatives(basis, dim)
-    return create_basis(name, description, X, basis, dbasis)
-end
diff --git a/src/lagrange_hexahedrons.jl b/src/lagrange_hexahedrons.jl
index 1e271ad..f152f8b 100644
--- a/src/lagrange_hexahedrons.jl
+++ b/src/lagrange_hexahedrons.jl
@@ -14,7 +14,7 @@ code = create_basis(
      ( 1.0,  1.0,  1.0), # N7
      (-1.0,  1.0,  1.0), # N8
     ),
-    "1 + u + v + w + u*v + v*w + w*u + u*v*w",
+    :(1 + u + v + w + u*v + v*w + w*u + u*v*w),
    )
 eval(code)
 
@@ -43,7 +43,7 @@ code = create_basis(
      ( 0.0,  1.0,  1.0), # N19
      (-1.0,  0.0,  1.0), # N20
     ),
-    "1 + u + v + w + u*v + v*w + u*w + u*v*w + u^2 + v^2 + w^2 + u^2*v + u*v^2 + v^2*w + v*w^2 + u*w^2 + u^2*w + u^2*v*w + u*v^2*w + u*v*w^2",
+    :(1 + u + v + w + u*v + v*w + u*w + u*v*w + u^2 + v^2 + w^2 + u^2*v + u*v^2 + v^2*w + v*w^2 + u*w^2 + u^2*w + u^2*v*w + u*v^2*w + u*v*w^2),
    )
 eval(code)
 
@@ -79,6 +79,6 @@ code = create_basis(
      ( 0.0,  0.0,  1.0), # N26
      ( 0.0,  0.0,  0.0), # N27
     ),
-    "1 + u + v + w + u*v + v*w + u*w + u*v*w + u^2 + v^2 + w^2 + u^2*v + u*v^2 + v^2*w + v*w^2 + u*w^2 + u^2*w + u^2*v*w + u*v^2*w + u*v*w^2 + u^2*v^2 + v^2*w^2 + u^2*w^2 + u^2*v^2*w + u*v^2*w^2 + u^2*v*w^2 + u^2*v^2*w^2",
+    :(1 + u + v + w + u*v + v*w + u*w + u*v*w + u^2 + v^2 + w^2 + u^2*v + u*v^2 + v^2*w + v*w^2 + u*w^2 + u^2*w + u^2*v*w + u*v^2*w + u*v*w^2 + u^2*v^2 + v^2*w^2 + u^2*w^2 + u^2*v^2*w + u*v^2*w^2 + u^2*v*w^2 + u^2*v^2*w^2),
    )
 eval(code)
diff --git a/src/lagrange_pyramids.jl b/src/lagrange_pyramids.jl
index 00d1bcc..73e0838 100644
--- a/src/lagrange_pyramids.jl
+++ b/src/lagrange_pyramids.jl
@@ -13,11 +13,11 @@ code = create_basis(
      ( 0.0,  0.0, 1.0), # N5
     ),
     (
-     "1/4*( (1+u)*(1+v) - w + u*v*w/(1-w) )",
-     "1/4*( (1+u)*(1-v) - w + u*v*w/(1-w) )",
-     "1/4*( (1-u)*(1-v) - w + u*v*w/(1-w) )",
-     "1/4*( (1-u)*(1+v) - w + u*v*w/(1-w) )",
-     "w",
+     :(1/4*( (1+u)*(1+v) - w + u*v*w/(1-w) )),
+     :(1/4*( (1+u)*(1-v) - w + u*v*w/(1-w) )),
+     :(1/4*( (1-u)*(1-v) - w + u*v*w/(1-w) )),
+     :(1/4*( (1-u)*(1+v) - w + u*v*w/(1-w) )),
+     :(1.0*w),
     ),
    )
 eval(code)
@@ -34,11 +34,11 @@ code = create_basis(
      ( 0.0,  0.0,  1.0), # N5
     ),
     (
-     "1/8 * (1-u) * (1-v) * (1-w)",
-     "1/8 * (1+u) * (1-v) * (1-w)",
-     "1/8 * (1+u) * (1+v) * (1-w)",
-     "1/8 * (1-u) * (1+v) * (1-w)",
-     "1/2 * (1+w)",
+     :(1/8 * (1-u) * (1-v) * (1-w)),
+     :(1/8 * (1+u) * (1-v) * (1-w)),
+     :(1/8 * (1+u) * (1+v) * (1-w)),
+     :(1/8 * (1-u) * (1+v) * (1-w)),
+     :(1/2 * (1+w)),
     ),
    )
 eval(code)
diff --git a/src/lagrange_quadrangles.jl b/src/lagrange_quadrangles.jl
index f7abfa0..ec37baa 100644
--- a/src/lagrange_quadrangles.jl
+++ b/src/lagrange_quadrangles.jl
@@ -10,7 +10,7 @@ code = create_basis(
      ( 1.0,  1.0), # N3
      (-1.0,  1.0)  # N4
     ),
-    "1 + u + v + u*v",
+    :(1 + u + v + u*v),
    )
 eval(code)
 
@@ -27,7 +27,7 @@ code = create_basis(
      ( 0.0,  1.0), # N7
      (-1.0,  0.0)  # N8
     ),
-    "1 + u + v + u*v + u^2 + u^2*v + u*v^2 + v^2",
+    :(1 + u + v + u*v + u^2 + u^2*v + u*v^2 + v^2),
    )
 eval(code)
 
@@ -45,6 +45,6 @@ code = create_basis(
      (-1.0,  0.0), # N8
      ( 0.0,  0.0)  # N9
     ),
-    "1 + u + v + u*v + u^2 + u^2*v + u*v^2 + v^2 + u^2*v^2",
+    :(1 + u + v + u*v + u^2 + u^2*v + u*v^2 + v^2 + u^2*v^2),
    )
 eval(code)
diff --git a/src/lagrange_segments.jl b/src/lagrange_segments.jl
index 3f43559..4bd93ac 100644
--- a/src/lagrange_segments.jl
+++ b/src/lagrange_segments.jl
@@ -5,12 +5,12 @@ code = create_basis(
     :Seg2,
     "2 node linear segment/line element",
     ( (-1.0,), ( 1.0,) ),
-    "1 + u")
+    :(1 + u))
 eval(code)
 
 code = create_basis(
     :Seg3,
     "3 node quadratic segment/line element",
     ( (-1.0,), (1.0,), (0.0,)),
-    "1 + u + u^2")
+    :(1 + u + u^2))
 eval(code)
diff --git a/src/lagrange_tetrahedrons.jl b/src/lagrange_tetrahedrons.jl
index 9f567ac..1a23047 100644
--- a/src/lagrange_tetrahedrons.jl
+++ b/src/lagrange_tetrahedrons.jl
@@ -10,7 +10,7 @@ code = create_basis(
      (0.0, 1.0, 0.0), # N3
      (0.0, 0.0, 1.0), # N4
     ),
-    "1 + u + v + w",
+    :(1 + u + v + w),
    )
 eval(code)
 
@@ -29,6 +29,6 @@ code = create_basis(
      (0.5, 0.0, 0.5), # N9
      (0.0, 0.5, 0.5), # N10
     ),
-    "1 + u + v + w + u*v + v*w + w*u + u^2 + v^2 + w^2",
+    :(1 + u + v + w + u*v + v*w + w*u + u^2 + v^2 + w^2),
    )
 eval(code)
diff --git a/src/lagrange_triangles.jl b/src/lagrange_triangles.jl
index eeee69f..c0cce0c 100644
--- a/src/lagrange_triangles.jl
+++ b/src/lagrange_triangles.jl
@@ -9,7 +9,7 @@ code = create_basis(
      (1.0, 0.0), # N2
      (0.0, 1.0), # N3
     ),
-    "1 + u + v",
+    :(1 + u + v),
    )
 eval(code)
 
@@ -24,7 +24,7 @@ code = create_basis(
      (0.5, 0.5), # N5
      (0.0, 0.5), # N6
     ),
-    "1 + u + v + u^2 + u*v + v^2",
+    :(1 + u + v + u^2 + u*v + v^2),
    )
 eval(code)
 
@@ -40,6 +40,6 @@ code = create_basis(
      (0.0, 0.5), # N6
      (1/3, 1/3), # N7
     ),
-    "1 + u + v + u^2 + u*v + v^2 + u^2*v^2",
+    :(1 + u + v + u^2 + u*v + v^2 + u^2*v^2),
    )
 eval(code)
diff --git a/src/lagrange_wedges.jl b/src/lagrange_wedges.jl
index d20a07b..7507ec1 100644
--- a/src/lagrange_wedges.jl
+++ b/src/lagrange_wedges.jl
@@ -14,12 +14,12 @@ code = create_basis(
      (0.0, 1.0,  1.0), # N6
     ),
     (
-     "1/2 * (1-w) * (1-u-v)", # N1
-     "1/2 * (1-w) * u", # N2
-     "1/2 * (1-w) * v", # N3
-     "1/2 * (1+w) * (1-u-v)", # N4
-     "1/2 * (1+w) * u", # N5
-     "1/2 * (1+w) * v", # N6
+     :(1/2 * (1-w) * (1-u-v)), # N1
+     :(1/2 * (1-w) * u), # N2
+     :(1/2 * (1-w) * v), # N3
+     :(1/2 * (1+w) * (1-u-v)), # N4
+     :(1/2 * (1+w) * u), # N5
+     :(1/2 * (1+w) * v), # N6
     ),
    )
 eval(code)
@@ -46,21 +46,21 @@ code = create_basis(
      (0.0, 1.0,  0.0), # N15
     ),
     (
-     "u^2*w^2 - u^2*w + 2*u*v*w^2 - 2*u*v*w - 3*u*w^2/2 + 3*u*w/2 + v^2*w^2 - v^2*w - 3*v*w^2/2 + 3*v*w/2 + w^2/2 - w/2",
-     "u^2*w^2 - u^2*w - u*w^2/2 + u*w/2",
-     "v^2*w^2 - v^2*w - v*w^2/2 + v*w/2",
-     "u^2*w^2 + u^2*w + 2*u*v*w^2 + 2*u*v*w - 3*u*w^2/2 - 3*u*w/2 + v^2*w^2 + v^2*w - 3*v*w^2/2 - 3*v*w/2 + w^2/2 + w/2",
-     "u^2*w^2 + u^2*w - u*w^2/2 - u*w/2",
-     "v^2*w^2 + v^2*w - v*w^2/2 - v*w/2",
-     "-2*u^2*w^2 + 2*u^2*w - 2*u*v*w^2 + 2*u*v*w + 2*u*w^2 - 2*u*w",
-     "2*u*v*w^2 - 2*u*v*w",
-     "-2*u*v*w^2 + 2*u*v*w - 2*v^2*w^2 + 2*v^2*w + 2*v*w^2 - 2*v*w",
-     "2*u^2*w^2 - 2*u^2*w - 4*u^2 + 2*u*v*w^2 - 2*u*v*w - 4*u*v - 2*u*w^2 + 2*u*w + 4*u",
-     "-2*u*v*w^2 + 2*u*v*w + 4*u*v",
-     "2*u*v*w^2 - 2*u*v*w - 4*u*v + 2*v^2*w^2 - 2*v^2*w - 4*v^2 - 2*v*w^2 + 2*v*w + 4*v",
-     "-2*u^2*w^2 + 2*u^2 - 4*u*v*w^2 + 4*u*v + 3*u*w^2 - 3*u - 2*v^2*w^2 + 2*v^2 + 3*v*w^2 - 3*v - w^2 + 1",
-     "-2*u^2*w^2 + 2*u^2 + u*w^2 - u",
-     "-2*v^2*w^2 + 2*v^2 + v*w^2 - v",
+     :(u^2*w^2 - u^2*w + 2*u*v*w^2 - 2*u*v*w - 3*u*w^2/2 + 3*u*w/2 + v^2*w^2 - v^2*w - 3*v*w^2/2 + 3*v*w/2 + w^2/2 - w/2),
+     :(u^2*w^2 - u^2*w - u*w^2/2 + u*w/2),
+     :(v^2*w^2 - v^2*w - v*w^2/2 + v*w/2),
+     :(u^2*w^2 + u^2*w + 2*u*v*w^2 + 2*u*v*w - 3*u*w^2/2 - 3*u*w/2 + v^2*w^2 + v^2*w - 3*v*w^2/2 - 3*v*w/2 + w^2/2 + w/2),
+     :(u^2*w^2 + u^2*w - u*w^2/2 - u*w/2),
+     :(v^2*w^2 + v^2*w - v*w^2/2 - v*w/2),
+     :(-2*u^2*w^2 + 2*u^2*w - 2*u*v*w^2 + 2*u*v*w + 2*u*w^2 - 2*u*w),
+     :(2*u*v*w^2 - 2*u*v*w),
+     :(-2*u*v*w^2 + 2*u*v*w - 2*v^2*w^2 + 2*v^2*w + 2*v*w^2 - 2*v*w),
+     :(2*u^2*w^2 - 2*u^2*w - 4*u^2 + 2*u*v*w^2 - 2*u*v*w - 4*u*v - 2*u*w^2 + 2*u*w + 4*u),
+     :(-2*u*v*w^2 + 2*u*v*w + 4*u*v),
+     :(2*u*v*w^2 - 2*u*v*w - 4*u*v + 2*v^2*w^2 - 2*v^2*w - 4*v^2 - 2*v*w^2 + 2*v*w + 4*v),
+     :(-2*u^2*w^2 + 2*u^2 - 4*u*v*w^2 + 4*u*v + 3*u*w^2 - 3*u - 2*v^2*w^2 + 2*v^2 + 3*v*w^2 - 3*v - w^2 + 1),
+     :(-2*u^2*w^2 + 2*u^2 + u*w^2 - u),
+     :(-2*v^2*w^2 + 2*v^2 + v*w^2 - v),
     ),
    )
 eval(code)
diff --git a/src/math.jl b/src/math.jl
index 30122e7..6ca6f1f 100644
--- a/src/math.jl
+++ b/src/math.jl
@@ -246,7 +246,7 @@ function eval_basis!(bi::BasisInfo{B}, X, xi) where B
         bi.invJ[7] = 1.0 / bi.detJ * (d*h - e*g)
         bi.invJ[8] = 1.0 / bi.detJ * (b*g - a*h)
         bi.invJ[9] = 1.0 / bi.detJ * (a*e - b*d)
-        A_mul_B!(bi.grad, bi.invJ, bi.dN)
+        mul!(bi.grad, bi.invJ, bi.dN)
     elseif dims == (2, 2)
         a, b, c, d = bi.J
         bi.detJ = a*d - b*c
@@ -254,7 +254,7 @@ function eval_basis!(bi::BasisInfo{B}, X, xi) where B
         bi.invJ[2] = 1.0 / bi.detJ * -b
         bi.invJ[3] = 1.0 / bi.detJ * -c
         bi.invJ[4] = 1.0 / bi.detJ * a
-        A_mul_B!(bi.grad, bi.invJ, bi.dN)
+        mul!(bi.grad, bi.invJ, bi.dN)
     elseif dims == (1, 1)
         bi.detJ = bi.J[1]
         bi.invJ[1] = 1.0 / bi.detJ
diff --git a/src/subs.jl b/src/subs.jl
new file mode 100644
index 0000000..20b5f28
--- /dev/null
+++ b/src/subs.jl
@@ -0,0 +1,56 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+
+using Calculus
+
+function subs(p::Number, ::Any)
+    return p
+end
+
+function subs(p::Symbol, data::Pair{Symbol, T}) where T
+    k, v = data
+    if p == k
+        return v
+    end
+    return p
+end
+
+function subs(p::Symbol, data::NTuple{N,Pair{Symbol, T}}) where {N, T}
+    for (k, v) in data
+        if p == k
+            return v
+        end
+    end
+    return p
+end
+
+function subs(p::Expr, d::Pair)
+    v = copy(p)
+    for j in 2:length(p.args)
+        v.args[j] = subs(v.args[j], d)
+    end
+    return v
+end
+
+"""
+    subs(expression, data)
+
+Given expression and pair(s) of `symbol => value` data, substitute to expression.
+
+# Examples
+
+Let us have polynomial `1 + u + v + u*v^2`, and substitute u=1 and v=2:
+```julia
+expression = :(1 + u + v + u*v^2)
+data = (:u => 1.0, :v => 2.0)
+subs(expression, data)
+8.0
+```
+
+"""
+function subs(p::Expr, data::NTuple{N,Pair{Symbol, T}}) where {N, T}
+    for di in data
+        p = subs(p, di)
+    end
+    return simplify(p)
+end
diff --git a/src/vandermonde.jl b/src/vandermonde.jl
new file mode 100644
index 0000000..16455e8
--- /dev/null
+++ b/src/vandermonde.jl
@@ -0,0 +1,52 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+
+"""
+    vandermonde_matrix(polynomial, coordinates)
+
+Given some polynomial and coordinates points (1-3 dimensions), create a Vandermonde
+matrix.
+
+# Example
+
+To genererate a Vandermonde matrix for a reference quadrangle `[-1.0, 1.0]^2` for
+polynomial `p(u,v) = 1 + u + v + u*v`, one writes:
+
+```julia
+polynomial = :(1 + u + v + u*v)
+coordinates = ((-1.0,-1.0), (1.0,-1.0), (1.0,1.0), (-1.0,1.0))
+V = vandermonde_matrix(polynomial, coordinates)
+
+# output
+
+[
+1.0 -1.0 -1.0  1.0
+1.0  1.0 -1.0 -1.0
+1.0  1.0  1.0  1.0
+1.0 -1.0  1.0 -1.0
+]
+
+```
+
+# References
+- Wikipedia contributors. (2018, August 1). Vandermonde matrix. In Wikipedia, The Free Encyclopedia. Retrieved 10:00, August 20, 2018, from https://en.wikipedia.org/w/index.php?title=Vandermonde_matrix&oldid=852930962
+"""
+function vandermonde_matrix(polynomial::Expr, coordinates::NTuple{N, NTuple{D, T}}) where {N, D, T<:Number}
+    A = zeros(N, N)
+    first(polynomial.args) == :+ || error("Use only summation between terms of polynomial")
+    args = polynomial.args[2:end]
+    for i in 1:N
+        X = coordinates[i]
+        if D == 1
+            data = (:u => X[1],)
+        elseif D == 2
+            data = (:u => X[1], :v => X[2])
+        elseif D == 3
+            data = (:u => X[1], :v => X[2], :w => X[3])
+        end
+        for (j, term) in enumerate(args)
+            A[i,j] = subs(term, data)
+        end
+    end
+    return A
+end
diff --git a/test/REQUIRE b/test/REQUIRE
index 0ec8121..e69de29 100644
--- a/test/REQUIRE
+++ b/test/REQUIRE
@@ -1 +0,0 @@
-TimerOutputs
diff --git a/test/runtests.jl b/test/runtests.jl
index 4ad1c60..9199fe4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,23 +1,18 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
 
-using Base.Test
-using TimerOutputs
-const to = TimerOutput()
-
-test_files = String[]
-push!(test_files, "test_create_basis.jl")
-push!(test_files, "test_lagrange_elements.jl")
-push!(test_files, "test_nurbs_elements.jl")
-push!(test_files, "test_math.jl")
+if VERSION < v"0.7.0"
+    using Base.Test
+else
+    using Test
+end
 
 @testset "FEMBasis.jl" begin
-    for fn in test_files
-        timeit(to, fn) do
-            include(fn)
-        end
-    end
+    @testset "test substitution" begin include("test_substitute.jl") end
+    @testset "create Vandermonde matrix" begin include("test_vandermonde.jl") end
+    @testset "create interpolation polynomials" begin include("test_generate_polynomials.jl") end
+    @testset "test_create_basis" begin include("test_create_basis.jl") end
+    @testset "test_lagrange_elements" begin include("test_lagrange_elements.jl") end
+    @testset "test_nurbs_elements" begin include("test_nurbs_elements.jl") end
+    @testset "test_math" begin include("test_math.jl") end
 end
-println()
-println("Test statistics:")
-println(to)
diff --git a/test/test_create_basis.jl b/test/test_create_basis.jl
index 2936bdb..ac1a341 100644
--- a/test/test_create_basis.jl
+++ b/test/test_create_basis.jl
@@ -1,110 +1,72 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
 
-using Base.Test
 using FEMBasis
-#using FEMBasis: calculate_basis_coefficients, calculate_interpolation_polynomials,
-#                calculate_interpolation_polynomial_derivatives, create_lagrange_basis
+using FEMBasis: create_basis
 
-@testset "FEMBasis.jl" begin
-
-@testset "Calculate interpolation polynomial matrix" begin
-    p = :(1 + u + v)
-    X = ((0.0, 0.0), (1.0, 0.0), (0.0, 1.0))
-    A = FEMBasis.calculate_basis_coefficients(p, X)
-    A_expected = [1.0 0.0 0.0; 1.0 1.0 0.0; 1.0 0.0 1.0]
-    @test isapprox(A, A_expected)
-end
-
-@testset "Calculate interpolation polynomials" begin
-    p = :(1 + u + v)
-    A = [1.0 0.0 0.0; 1.0 1.0 0.0; 1.0 0.0 1.0]
-    equations = FEMBasis.calculate_interpolation_polynomials(p, A)
-    @test equations[1] == :(1.0 + -u + -v)
-    @test equations[2] == :(+u)
-    @test equations[3] == :(+v)
+if VERSION < v"0.7.0"
+    using Base.Test
+else
+    using Test
 end
 
+global X = ((0.0, 0.0), (1.0, 0.0), (0.0, 1.0))
+global basis = [:(1.0 - u - v), :(1.0u), :(1.0v)]
+global dbasis = Vector[[-1.0, -1.0], [1.0, 0.0], [0.0, 1.0]]
 
-@testset "Calculate derivatives of interpolation polynomials" begin
-    basis = [:(1.0 - u - v), :(1.0u), :(1.0v)]
-    dbasis = FEMBasis.calculate_interpolation_polynomial_derivatives(basis, 2)
-    @test isapprox(dbasis[1][1], -1.0)
-    @test isapprox(dbasis[1][2], -1.0)
-    @test isapprox(dbasis[2][1], 1.0)
-    @test isapprox(dbasis[2][2], 0.0)
-    @test isapprox(dbasis[3][1], 0.0)
-    @test isapprox(dbasis[3][2], 1.0)
+code = create_basis(:TestTriangle1, "test triangle 1", X, basis, dbasis)
+eval(code)
+N = zeros(1,size(TestTriangle1, 2))
+X = get_reference_element_coordinates(TestTriangle1)
+for i=1:length(TestTriangle1)
+    eval_basis!(TestTriangle1, N, X[i])
+    N_expected = zeros(1, length(TestTriangle1))
+    N_expected[i] = 1.0
+    @test isapprox(N, N_expected)
 end
+dN = zeros(size(TestTriangle1)...)
+eval_dbasis!(TestTriangle1, dN, X[1])
+@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
+eval_dbasis!(TestTriangle1, dN, X[2])
+@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
+eval_dbasis!(TestTriangle1, dN, X[3])
+@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
 
-@testset "test create basis, N and dN given" begin
-    X = ((0.0, 0.0), (1.0, 0.0), (0.0, 1.0))
-    basis = [:(1.0 - u - v), :(1.0u), :(1.0v)]
-    dbasis = Vector[[-1.0, -1.0], [1.0, 0.0], [0.0, 1.0]]
-    code = FEMBasis.create_basis(:TestTriangle, "test triangle", X, basis, dbasis)
-    eval(code)
-    N = zeros(1,size(TestTriangle, 2))
-    X = get_reference_element_coordinates(TestTriangle)
-    for i=1:length(TestTriangle)
-        eval_basis!(TestTriangle, N, X[i])
-        N_expected = zeros(1, length(TestTriangle))
-        N_expected[i] = 1.0
-        @test isapprox(N, N_expected)
-    end
-    dN = zeros(size(TestTriangle)...)
-    eval_dbasis!(TestTriangle, dN, X[1])
-    @test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
-    eval_dbasis!(TestTriangle, dN, X[2])
-    @test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
-    eval_dbasis!(TestTriangle, dN, X[3])
-    @test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
-end
 
-@testset "test create basis, N given" begin
-    X = ((0.0, 0.0), (1.0, 0.0), (0.0, 1.0))
-    basis = ( :(1 - u - v), :(1.0u), :(1.0v) )
-    code = FEMBasis.create_basis(:TestTriangle, "test triangle", X, basis)
-    eval(code)
-    N = zeros(1,size(TestTriangle, 2))
-    X = get_reference_element_coordinates(TestTriangle)
-    for i=1:length(TestTriangle)
-        eval_basis!(TestTriangle, N, X[i])
-        N_expected = zeros(1, length(TestTriangle))
-        N_expected[i] = 1.0
-        @test isapprox(N, N_expected)
-    end
-    dN = zeros(size(TestTriangle)...)
-    eval_dbasis!(TestTriangle, dN, X[1])
-    @test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
-    eval_dbasis!(TestTriangle, dN, X[2])
-    @test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
-    eval_dbasis!(TestTriangle, dN, X[3])
-    @test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
+code = create_basis(:TestTriangle2, "test triangle 2", X, basis)
+eval(code)
+N = zeros(1,size(TestTriangle2, 2))
+X = get_reference_element_coordinates(TestTriangle2)
+for i=1:length(TestTriangle2)
+    eval_basis!(TestTriangle2, N, X[i])
+    N_expected = zeros(1, length(TestTriangle2))
+    N_expected[i] = 1.0
+    @test isapprox(N, N_expected)
 end
+dN = zeros(size(TestTriangle2)...)
+eval_dbasis!(TestTriangle2, dN, X[1])
+@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
+eval_dbasis!(TestTriangle2, dN, X[2])
+@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
+eval_dbasis!(TestTriangle2, dN, X[3])
+@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
 
-@testset "test create basis, ansatz polynomial given" begin
-    X = ((0.0, 0.0), (1.0, 0.0), (0.0, 1.0))
-    p1 = :(1 + u + v)
-    p2 = "1 + u + v"
-    code1 = FEMBasis.create_basis(:TestTriangle, "test triangle", X, p1)
-    code2 = FEMBasis.create_basis(:TestTriangle, "test triangle", X, p2)
-    @test code1 == code2
-    eval(code1)
-    N = zeros(1,size(TestTriangle, 2))
-    X = get_reference_element_coordinates(TestTriangle)
-    for i=1:length(TestTriangle)
-        eval_basis!(TestTriangle, N, X[i])
-        N_expected = zeros(1, length(TestTriangle))
-        N_expected[i] = 1.0
-        @test isapprox(N, N_expected)
-    end
-    dN = zeros(size(TestTriangle)...)
-    eval_dbasis!(TestTriangle, dN, X[1])
-    @test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
-    eval_dbasis!(TestTriangle, dN, X[2])
-    @test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
-    eval_dbasis!(TestTriangle, dN, X[3])
-    @test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
-end
 
+p = :(1 + u + v)
+code = create_basis(:TestTriangle3, "test triangle 3", X, p)
+eval(code)
+N = zeros(1, size(TestTriangle3, 2))
+X = get_reference_element_coordinates(TestTriangle3)
+for i=1:length(TestTriangle3)
+    eval_basis!(TestTriangle3, N, X[i])
+    N_expected = zeros(1, length(TestTriangle3))
+    N_expected[i] = 1.0
+    @test isapprox(N, N_expected)
 end
+dN = zeros(size(TestTriangle3)...)
+eval_dbasis!(TestTriangle3, dN, X[1])
+@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
+eval_dbasis!(TestTriangle3, dN, X[2])
+@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
+eval_dbasis!(TestTriangle3, dN, X[3])
+@test isapprox(dN, [-1.0 1.0 0.0; -1.0 0.0 1.0])
diff --git a/test/test_generate_polynomials.jl b/test/test_generate_polynomials.jl
new file mode 100644
index 0000000..05f5442
--- /dev/null
+++ b/test/test_generate_polynomials.jl
@@ -0,0 +1,20 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+
+using FEMBasis: calculate_interpolation_polynomials,
+                calculate_interpolation_polynomial_derivatives
+
+p = :(1 + u + v)
+A = [1.0 0.0 0.0; 1.0 1.0 0.0; 1.0 0.0 1.0]
+basis = calculate_interpolation_polynomials(p, A)
+@test basis[1] == :(1.0 + -u + -v)
+@test basis[2] == :(+u)
+@test basis[3] == :(+v)
+
+dbasis = calculate_interpolation_polynomial_derivatives(basis, 2)
+@test isapprox(dbasis[1][1], -1.0)
+@test isapprox(dbasis[1][2], -1.0)
+@test isapprox(dbasis[2][1], 1.0)
+@test isapprox(dbasis[2][2], 0.0)
+@test isapprox(dbasis[3][1], 0.0)
+@test isapprox(dbasis[3][2], 1.0)
diff --git a/test/test_lagrange_elements.jl b/test/test_lagrange_elements.jl
index e5fca59..1eb14d9 100644
--- a/test/test_lagrange_elements.jl
+++ b/test/test_lagrange_elements.jl
@@ -1,7 +1,11 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
 
-using Base.Test
+if VERSION < v"0.7.0"
+    using Base.Test
+else
+    using Test
+end
 
 using FEMBasis
 
diff --git a/test/test_math.jl b/test/test_math.jl
index 5e149e6..715c1e0 100644
--- a/test/test_math.jl
+++ b/test/test_math.jl
@@ -1,125 +1,122 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
 
-using Base.Test
+if VERSION < v"0.7.0"
+    using Base.Test
+else
+    using Test
+end
+
 using FEMBasis
+using LinearAlgebra
 
-@testset "interpolate scalar field" begin
-    # in unit square: T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
-    B = Quad4()
-    X = ((0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0))
-    T = (1.0, 2.0, 3.0, 4.0)
-    T_known(X) = 1 + X[1] + 3*X[2] - 2*X[1]*X[2]
-    T_interpolated = interpolate(B, T, (0.0, 0.0))
-    @test isapprox(T_interpolated, T_known((0.5, 0.5)))
-end
+# interpolate scalar field in unit square:
+# T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
+B = Quad4()
+X = ((0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0))
+T = (1.0, 2.0, 3.0, 4.0)
+T_known(X) = 1 + X[1] + 3*X[2] - 2*X[1]*X[2]
+T_interpolated = interpolate(B, T, (0.0, 0.0))
+@test isapprox(T_interpolated, T_known((0.5, 0.5)))
 
-@testset "calculate gradient dB/dX" begin
-    B = Quad4()
-    X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
-    dBdX = grad(B, X, (0.0, 0.0))
-    @test isapprox(dBdX, 1/2*[-1 1 1 -1; -1 -1 1 1])
-end
+# calculate gradient dB/dX
+B = Quad4()
+X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
+dBdX = grad(B, X, (0.0, 0.0))
+@test isapprox(dBdX, 1/2*[-1 1 1 -1; -1 -1 1 1])
 
-@testset "interpolate gradient of scalar field" begin
-    # in unit square: grad(T)(X) = [1-2X[2], 3-2*X[1]]
-    B = Quad4()
-    X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
-    T = (1.0, 2.0, 3.0, 4.0)
-    dTdX = grad(B, T, X, (0.0, 0.0))
-    dTdX_expected(X) = [1-2*X[2] 3-2*X[1]]
-    @test isapprox(dTdX, dTdX_expected([0.5, 0.5]))
-end
+# interpolate gradient of scalar field in unit square:
+# grad(T)(X) = [1-2X[2], 3-2*X[1]]
+B = Quad4()
+X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
+T = (1.0, 2.0, 3.0, 4.0)
+dTdX = grad(B, T, X, (0.0, 0.0))
+dTdX_expected(X) = [1-2*X[2] 3-2*X[1]]
+@test isapprox(dTdX, dTdX_expected([0.5, 0.5]))
 
-@testset "interpolate vector field" begin
-    # in unit square, u(X,t) = [1/4*t*X[1]*X[2], 0, 0]
-    B = Quad4()
-    u = Vector[[0.0, 0.0], [0.0, 0.0], [1/4, 0.0], [0.0, 0.0]]
-    u_known(X) = [1/4*X[1]*X[2], 0]
-    u_interpolated = interpolate(B, u, (0.0, 0.0))
-    @test isapprox(u_interpolated, u_known((0.5, 0.5)))
-end
+# interpolate vector field in unit square:
+# u(X,t) = [1/4*t*X[1]*X[2], 0, 0]
+B = Quad4()
+u = Vector[[0.0, 0.0], [0.0, 0.0], [1/4, 0.0], [0.0, 0.0]]
+u_known(X) = [1/4*X[1]*X[2], 0]
+u_interpolated = interpolate(B, u, (0.0, 0.0))
+@test isapprox(u_interpolated, u_known((0.5, 0.5)))
 
-@testset "interpolate gradient of vector field" begin
-    # in unit square, u(X) = t*[X[1]*(X[2]+1), X[1]*(4*X[2]-1)]
-    # => u_i,j = t*[X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
-    B = Quad4()
-    X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
-    u = ([0.0, 0.0], [1.0, -1.0], [2.0, 3.0], [0.0, 0.0])
-    dudX = grad(B, u, X, (0.0, 0.0))
-    dudX_expected(X) = [X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
-    @test isapprox(dudX, dudX_expected([0.5, 0.5]))
-end
+# interpolate gradient of vector field in unit square:
+# u(X) = t*[X[1]*(X[2]+1), X[1]*(4*X[2]-1)]
+# => u_i,j = t*[X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
+B = Quad4()
+X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
+u = ([0.0, 0.0], [1.0, -1.0], [2.0, 3.0], [0.0, 0.0])
+dudX = grad(B, u, X, (0.0, 0.0))
+dudX_expected(X) = [X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
+@test isapprox(dudX, dudX_expected([0.5, 0.5]))
 
-@testset "test BasisInfo" begin
-    B = BasisInfo(Seg2)
-    @test length(B) == 2
-    @test size(B) == (1, 2)
-    X = ((0.0,), (1.1))
-    xi = (0.0,)
-    eval_basis!(B, X, xi)
-    @test isapprox(B.invJ, inv(B.J))
-    @test isapprox(B.detJ, det(B.J))
-    B = BasisInfo(Quad4)
-    X = ((0.0,0.0), (1.1,0.0), (1.0,1.0), (0.0,1.0))
-    xi = (0.0,0.0)
-    eval_basis!(B, X, xi)
-    @test isapprox(B.invJ, inv(B.J))
-    @test isapprox(B.detJ, det(B.J))
-    B = BasisInfo(Hex8)
-    X = ((0.0,0.0,0.0), (1.1,0.0,0.0), (1.0,1.0,0.0), (0.0,1.0,0.0),
-         (0.0,0.0,1.0), (1.0,0.0,1.0), (1.0,1.0,1.0), (0.0,1.0,1.0))
-    xi = (0.0,0.0,0.0)
-    eval_basis!(B, X, xi)
-    @test isapprox(B.invJ, inv(B.J))
-    @test isapprox(B.detJ, det(B.J))
-end
+# test BasisInfo
+B = BasisInfo(Seg2)
+@test length(B) == 2
+@test size(B) == (1, 2)
+X = ((0.0,), (1.1))
+xi = (0.0,)
+eval_basis!(B, X, xi)
+@test isapprox(B.invJ, inv(B.J))
+@test isapprox(B.detJ, det(B.J))
+B = BasisInfo(Quad4)
+X = ((0.0,0.0), (1.1,0.0), (1.0,1.0), (0.0,1.0))
+xi = (0.0,0.0)
+eval_basis!(B, X, xi)
+@test isapprox(B.invJ, inv(B.J))
+@test isapprox(B.detJ, det(B.J))
+B = BasisInfo(Hex8)
+X = ((0.0,0.0,0.0), (1.1,0.0,0.0), (1.0,1.0,0.0), (0.0,1.0,0.0),
+     (0.0,0.0,1.0), (1.0,0.0,1.0), (1.0,1.0,1.0), (0.0,1.0,1.0))
+xi = (0.0,0.0,0.0)
+eval_basis!(B, X, xi)
+@test isapprox(B.invJ, inv(B.J))
+@test isapprox(B.detJ, det(B.J))
 
-@testset "evaluate gradient using BasisInfo" begin
-    B = BasisInfo(Quad4)
-    X = ((0.0,0.0), (1.0,0.0), (1.0,1.0), (0.0,1.0))
-    xi = (0.0, 0.0)
-    eval_basis!(B, X, xi)
-    u = ((0.0, 0.0), (1.0, -1.0), (2.0, 3.0), (0.0, 0.0))
-    gradu = zeros(2, 2)
-    grad!(B, gradu, u)
-    gradu_expected = [1.5 0.5; 1.0  2.0]
-    @test isapprox(gradu, gradu_expected)
-end
+# evaluate gradient using BasisInfo
+B = BasisInfo(Quad4)
+X = ((0.0,0.0), (1.0,0.0), (1.0,1.0), (0.0,1.0))
+xi = (0.0, 0.0)
+eval_basis!(B, X, xi)
+u = ((0.0, 0.0), (1.0, -1.0), (2.0, 3.0), (0.0, 0.0))
+gradu = zeros(2, 2)
+grad!(B, gradu, u)
+gradu_expected = [1.5 0.5; 1.0  2.0]
+@test isapprox(gradu, gradu_expected)
 
-@testset "test curves" begin
-    bi = BasisInfo(Seg2)
-    X1 = ((0.0,0.0,0.0), (1.0,0.0,0.0))
-    X2 = ((0.0,0.0,0.0), (0.0,1.0,0.0))
-    X3 = ((0.0,0.0,0.0), (0.0,0.0,1.0))
-    xi = (0.0, 0.0)
-    eval_basis!(bi, X1, xi)
-    @test isapprox(bi.detJ, 0.5)
-    eval_basis!(bi, X2, xi)
-    @test isapprox(bi.detJ, 0.5)
-    eval_basis!(bi, X3, xi)
-    @test isapprox(bi.detJ, 0.5)
-    @test isapprox(jacobian(Seg2(), X1, xi), [0.5 0.0 0.0])
-    @test isapprox(jacobian(Seg2(), X2, xi), [0.0 0.5 0.0])
-    @test isapprox(jacobian(Seg2(), X3, xi), [0.0 0.0 0.5])
-end
+# test curves
+bi = BasisInfo(Seg2)
+X1 = ((0.0,0.0,0.0), (1.0,0.0,0.0))
+X2 = ((0.0,0.0,0.0), (0.0,1.0,0.0))
+X3 = ((0.0,0.0,0.0), (0.0,0.0,1.0))
+xi = (0.0, 0.0)
+eval_basis!(bi, X1, xi)
+@test isapprox(bi.detJ, 0.5)
+eval_basis!(bi, X2, xi)
+@test isapprox(bi.detJ, 0.5)
+eval_basis!(bi, X3, xi)
+@test isapprox(bi.detJ, 0.5)
+@test isapprox(jacobian(Seg2(), X1, xi), [0.5 0.0 0.0])
+@test isapprox(jacobian(Seg2(), X2, xi), [0.0 0.5 0.0])
+@test isapprox(jacobian(Seg2(), X3, xi), [0.0 0.0 0.5])
 
-@testset "test manifolds" begin
-    bi = BasisInfo(Quad4)
-    X1 = ((0.0,0.0,0.0), (1.0,0.0,0.0), (1.0,1.0,0.0), (0.0,1.0,0.0))
-    X2 = ((0.0,0.0,0.0), (0.0,1.0,0.0), (0.0,1.0,1.0), (0.0,0.0,1.0))
-    X3 = ((0.0,0.0,0.0), (0.0,0.0,1.0), (1.0,0.0,1.0), (1.0,0.0,0.0))
-    xi = (0.0, 0.0)
-    eval_basis!(bi, X1, xi)
-    @test isapprox(bi.detJ, 0.25)
-    eval_basis!(bi, X2, xi)
-    @test isapprox(bi.detJ, 0.25)
-    eval_basis!(bi, X3, xi)
-    @test isapprox(bi.detJ, 0.25)
-    J1 = jacobian(Quad4(), X1, xi)
-    J2 = jacobian(Quad4(), X2, xi)
-    J3 = jacobian(Quad4(), X3, xi)
-    @test isapprox(J1, [0.5 0.0 0.0; 0.0 0.5 0.0])
-    @test isapprox(J2, [0.0 0.5 0.0; 0.0 0.0 0.5])
-    @test isapprox(J3, [0.0 0.0 0.5; 0.5 0.0 0.0])
-end
+# test manifolds
+bi = BasisInfo(Quad4)
+X1 = ((0.0,0.0,0.0), (1.0,0.0,0.0), (1.0,1.0,0.0), (0.0,1.0,0.0))
+X2 = ((0.0,0.0,0.0), (0.0,1.0,0.0), (0.0,1.0,1.0), (0.0,0.0,1.0))
+X3 = ((0.0,0.0,0.0), (0.0,0.0,1.0), (1.0,0.0,1.0), (1.0,0.0,0.0))
+xi = (0.0, 0.0)
+eval_basis!(bi, X1, xi)
+@test isapprox(bi.detJ, 0.25)
+eval_basis!(bi, X2, xi)
+@test isapprox(bi.detJ, 0.25)
+eval_basis!(bi, X3, xi)
+@test isapprox(bi.detJ, 0.25)
+J1 = jacobian(Quad4(), X1, xi)
+J2 = jacobian(Quad4(), X2, xi)
+J3 = jacobian(Quad4(), X3, xi)
+@test isapprox(J1, [0.5 0.0 0.0; 0.0 0.5 0.0])
+@test isapprox(J2, [0.0 0.5 0.0; 0.0 0.0 0.5])
+@test isapprox(J3, [0.0 0.0 0.5; 0.5 0.0 0.0])
diff --git a/test/test_nurbs_elements.jl b/test/test_nurbs_elements.jl
index 896de37..bf1e300 100644
--- a/test/test_nurbs_elements.jl
+++ b/test/test_nurbs_elements.jl
@@ -1,32 +1,31 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
 
-using Base.Test
+if VERSION < v"0.7.0"
+    using Base.Test
+else
+    using Test
+end
+
 using FEMBasis
 
-@testset "NSeg interpolate" begin
-    basis = NSeg()
-    @test size(basis) == (1, 2)
-    @test length(basis) == 2
-    N = zeros(1, 2)
-    eval_basis!(basis, N, (0.0, ))
-    @test isapprox(N, [0.5 0.5])
-end
+basis = NSeg()
+@test size(basis) == (1, 2)
+@test length(basis) == 2
+N = zeros(1, 2)
+eval_basis!(basis, N, (0.0, ))
+@test isapprox(N, [0.5 0.5])
 
-@testset "NSurf interpolate" begin
-    basis = NSurf()
-    @test size(basis) == (2, 4)
-    @test length(basis) == 4
-    N = zeros(1, 4)
-    eval_basis!(basis, N, (0.0, 0.0))
-    @test isapprox(N, [0.25 0.25 0.25 0.25])
-end
+basis = NSurf()
+@test size(basis) == (2, 4)
+@test length(basis) == 4
+N = zeros(1, 4)
+eval_basis!(basis, N, (0.0, 0.0))
+@test isapprox(N, [0.25 0.25 0.25 0.25])
 
-@testset "NSolid interpolate" begin
-    basis = NSolid()
-    @test size(basis) == (3, 8)
-    @test length(basis) == 8
-    N = zeros(1, 8)
-    eval_basis!(basis, N, (0.0, 0.0, 0.0))
-    @test isapprox(N, [0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125])
-end
+basis = NSolid()
+@test size(basis) == (3, 8)
+@test length(basis) == 8
+N = zeros(1, 8)
+eval_basis!(basis, N, (0.0, 0.0, 0.0))
+@test isapprox(N, [0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125])
diff --git a/test/test_substitute.jl b/test/test_substitute.jl
new file mode 100644
index 0000000..83e4f5b
--- /dev/null
+++ b/test/test_substitute.jl
@@ -0,0 +1,15 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+
+using FEMBasis: subs
+if VERSION < v"0.7.0"
+    using Base.Test
+else
+    using Test
+end
+
+expression = :(1 + u + v + u*v^2)
+data = (:u => 1.0, :v => 2.0)
+result = subs(expression, data)
+@debug "subs(expression, data) where" expression data
+@test isapprox(result, 8.0)
diff --git a/test/test_vandermonde.jl b/test/test_vandermonde.jl
new file mode 100644
index 0000000..29247ab
--- /dev/null
+++ b/test/test_vandermonde.jl
@@ -0,0 +1,30 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+
+if VERSION < v"0.7.0"
+    using Base.Test
+else
+    using Test
+end
+
+using FEMBasis: vandermonde_matrix
+
+polynomial = :(1 + u + v + u*v)
+coordinates = ((-1.0,-1.0), (1.0,-1.0), (1.0,1.0), (-1.0,1.0))
+V = vandermonde_matrix(polynomial, coordinates)
+@debug "create Vandermonde matrix" polynomial coordinates V
+V_expected = [
+    1.0 -1.0 -1.0  1.0
+    1.0  1.0 -1.0 -1.0
+    1.0  1.0  1.0  1.0
+    1.0 -1.0  1.0 -1.0]
+@test isapprox(V, V_expected)
+    
+polynomial = :(1 + u + v)
+coordinates = ((0.0, 0.0), (1.0, 0.0), (0.0, 1.0))
+V = vandermonde_matrix(polynomial, coordinates)
+V_expected = [
+    1.0 0.0 0.0
+    1.0 1.0 0.0
+    1.0 0.0 1.0]
+@test isapprox(V, V_expected)