diff -r e958ed76e484 -r 0090a86d8b72 mononomial.m
--- a/mononomial.m	Tue Jun 26 16:59:03 2018 +0200
+++ b/mononomial.m	Wed Jun 27 11:10:52 2018 +0200
@@ -1,8 +1,17 @@
-function y = mononomial(x, k)
-    if k < 0
-        y = x*0;
+% calculate a N-D mononomial with powers k in points x:
+%  z = x(:,1).^k(1) * x(:,2).^k(2) * ...
+function z = mononomial(x, k)
+    assert(size(x,2) == length(k), 'k must have the same length as the width of x');
+
+    if any(k < 0)
+        z = x(:,1)*0;
         return
     end
-    y = x.^k/factorial(k);
+
+    denom = prod(factorial(k));
+
+    for i = 1:length(k)
+        x(:,i) = x(:,i).^k(i);
+    end
+    z = prod(x,2)/denom;
 end
-

diff -r e958ed76e484 -r 0090a86d8b72 vandermonde.m
--- a/vandermonde.m	Tue Jun 26 16:59:03 2018 +0200
+++ b/vandermonde.m	Wed Jun 27 11:10:52 2018 +0200
@@ -1,10 +1,15 @@
 % Create vandermonde matrix for points x and polynomials of order p
-% x and p are vectors
-% v is a length(x) by length(p) matrix
+% x is a list of N points of size [N,dim],
+% p is a list of polynomial orders of size [M, dim].
+% the given mononomials are evaluated and the NxM matrix V is returned.
 function V = vandermonde(x, p)
-    V = sym(zeros(length(x), length(p))); % Is there a way to make this work for both double and sym
+    assert(size(x,2) == size(p,2), 'x and p must have the same number of columns')
+    n = size(x,1);
+    m = size(p,1);
 
-    for i = 1:length(p)
-        V(:, i) = mononomial(x,p(i));
+    for i = 1:m
+        V(:,i) = mononomial(x, p(i,:));
     end
+
+    assertSize(V,[n,m]);
 end