--- a/mapFeature.m
+++ b/mapFeature.m
@@ -1,21 +1,21 @@
 function out = mapFeature(X1, X2)
-% MAPFEATURE Feature mapping function to polynomial features
-%
-%   MAPFEATURE(X1, X2) maps the two input features
-%   to quadratic features used in the regularization exercise.
-%
-%   Returns a new feature array with more features, comprising of 
-%   X1, X2, X1.^2, X2.^2, X1*X2, X1*X2.^2, etc..
-%
-%   Inputs X1, X2 must be the same size
-%
+  ## MAPFEATURE Feature mapping function to polynomial features
+  ##
+  ##   MAPFEATURE(X1, X2) maps the two input features
+  ##   to quadratic features used in the regularization exercise.
+  ##
+  ##   Returns a new feature array with more features, comprising of 
+  ##   X1, X2, X1.^2, X2.^2, X1*X2, X1*X2.^2, etc..
+  ##
+  ##   Inputs X1, X2 must be the same size
+  ##
 
-degree = 6;
-out = ones(size(X1(:,1)));
-for i = 1:degree
-    for j = 0:i
-        out(:, end+1) = (X1.^(i-j)).*(X2.^j);
-    end
-end
+  degree = 6;
 
-end
\ No newline at end of file
+  ## Compute the powers with an upper-triangular matrix trick
+  [i, j] = find (triu (ones (degree+1, degree+1)));
+  i--; j--;
+  [j, i] = deal (i', j');
+
+  out = bsxfun (@power, X1, i - j).*bsxfun (@power, X2, j);
+endfunction