diff --git a/pdeutils/functions/interpolateFunction.m b/pdeutils/functions/interpolateFunction.m
index c1f8a6a60689e99ced13b3866280fa00ffce2fa0..5162d477a0d375817f6725e9d4f0076d2bdf0b2a 100644
--- a/pdeutils/functions/interpolateFunction.m
+++ b/pdeutils/functions/interpolateFunction.m
@@ -1,30 +1,36 @@
 function u = interpolateFunction(grid, basis, varargin)
+%interpolateFunction Interpolates a function with respect to a given finite
+%element space.
+%
+% Input:
+%   grid        The grid of the finite element space.
+%   basis       A basis of the finite element space.
+%   [varargin]  A list of function handles that are required by the
+%               localInterpolation method of the supplied basis.
+%
+% Output:
+%   u           Vector representation of the finite element interpolation.
+
+    % Dimension of the ansatz space.
+    n   = basis.size(grid);
+
+    % Initialize vector.
+    u = zeros(n, 1);
+
+    % Loop over all elements.
+    for e = 1:size(grid.elements,2)
+
+        % Compute local interpolation.
+        uLocal = basis.localInterpolation(grid, e, varargin{:});
+
+        % Local size.
+        nLocal = basis.localSize(grid, e);
+
+        % Compute global indices.
+        for j=1:nLocal
+            u(basis.index(grid, e, j)) = uLocal(j);
+        end
 
-% dimension of ansatz space
-n = basis.size(grid);
-
-% initialize vector
-u = zeros(n, 1);
-
-% loop over all elements
-for e=1:size(grid.elements,2)
-
-    E = grid.nodes(1, grid.elements(:,e));
-
-    %  transform function and derivatives to element to 
-    fLocal = {};
-    for i=1:length(varargin)
-        fLocal{i} = @(x)varargin{i}(E(1) + (E(2)-E(1))*x);
-    end
-
-    % compute local interpolation
-    uLocal = basis.localInterpolation(grid, e, fLocal{:});
-
-    % local size
-    nLocal = basis.localSize(grid, e);
-
-    % compute global indices
-    for j=1:nLocal
-        u(basis.index(grid, e, j)) = uLocal(j);
     end
+    
 end