--- 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
-