--- a/src/Makefile.in	Sat Mar 01 14:11:24 2014 +0100
+++ b/src/Makefile.in	Tue Mar 04 12:19:04 2014 +0000
@@ -73,10 +73,10 @@
 	$(MKOCTFILE) $(CPPFLAGS) -c save.cc $(LDFLAGS) -o $@ -I.
 
 mkfunction:
-	 mkdir @function
+	 mkdir -p @function
 
 mkmesh:
-	 mkdir @mesh
+	 mkdir -p @mesh
 
 plot_mesh.oct: plot_mesh.o mkmesh
 	$(MKOCTFILE) $(CPPFLAGS) plot_mesh.o -o ./@mesh/plot.oct $(LDFLAGS) $(LIBS)

--- a/src/feval.cc	Sat Mar 01 14:11:24 2014 +0100
+++ b/src/feval.cc	Tue Mar 04 12:19:04 2014 +0000
@@ -1,69 +1,140 @@
-/*
- Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net>
-
- This program is free software; you can redistribute it and/or modify it under
- the terms of the GNU General Public License as published by the Free Software
- Foundation; either version 3 of the License, or (at your option) any later
- version.
-
- This program is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
- details.
-
- You should have received a copy of the GNU General Public License along with
- this program; if not, see <http://www.gnu.org/licenses/>.
-*/
-
-#include "function.h"
-
-DEFUN_DLD (feval, args, , "-*- texinfo -*-\n\
-@deftypefn {Function File} {[@var{value}]} = \
-feval (@var{function_name}, @var{Coordinate})\n\
-Evaluate a function at a specific point of the domain and return the value. \n\
-The input parameters are the function and the point where it has to\
-be evaluated.\n\
-@seealso{Function}\n\
-@end deftypefn")
-{
-
-  int nargin = args.length ();
-  octave_value retval=0;
-  
-  if (nargin < 2 || nargin > 2)
-    print_usage ();
-  else
-    {
-      if (! function_type_loaded)
-        {
-          function::register_type ();
-          function_type_loaded = true;
-          mlock ();
-        }
-
-      if (args(0).type_id () == function::static_type_id ())
-        {
-          const function & fspo =
-            static_cast<const function&> (args(0).get_rep ());
-          Array<double> point= args(1).array_value ();
-
-          if (!error_state)
-            {
-              const boost::shared_ptr<const dolfin::Function> 
-                & f = fspo.get_pfun ();
-              dim_vector dims;
-              dims.resize (2);
-              dims(0) = 1;
-              dims(1) = f->value_dimension (0);
-              Array<double> res (dims);
-              dolfin::Array<double> x(point.length (), point.fortran_vec ());
-              dolfin::Array<double> values(res.length (), res.fortran_vec ());
-
-              f->eval (values, x);
-
-              retval = octave_value (res);
-            }
-        }
-    }
-  return retval;
-}
+/*
+ Copyright (C) 2013-2014 Marco Vassallo <gedeone-octave@users.sourceforge.net>
+
+ This program is free software; you can redistribute it and/or modify it under
+ the terms of the GNU General Public License as published by the Free Software
+ Foundation; either version 3 of the License, or (at your option) any later
+ version.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ details.
+
+ You should have received a copy of the GNU General Public License along with
+ this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "function.h"
+
+/**
+ * Internal helper routine to evaluate the given function at a single
+ * point (as Array<double>).  Returned is the function value also as
+ * Array<double>, which gets stored into the given array.
+ * too, so that they don't have to be constructed.
+ * @param f Dolfin function handle to evaluate.
+ * @param pt Point where we want to evaluate f as Octave array of coordinates.
+ * @param res Store the value here, should already be of the correct size.
+ */
+static void
+evaluate (const dolfin::Function& f, const Array<double>& pt,
+          Array<double>& res)
+{
+  Array<double>& ptNonconst = const_cast<Array<double>&> (pt);
+  dolfin::Array<double> x(ptNonconst.length (), ptNonconst.fortran_vec ());
+  dolfin::Array<double> values(res.length (), res.fortran_vec ());
+
+  f.eval (values, x);
+}
+
+DEFUN_DLD (feval, args, , "-*- texinfo -*-\n\
+@deftypefn {Function File} {[@var{value}]} = \
+feval (@var{function_name}, @var{Coordinate})\n\
+@deftypefnx {Function File} {[@var{value}]} = \
+feval (@var{function_name}, @var{x}, @var{y})\n\
+@deftypefnx {Function File} {[@var{value}]} = \
+feval (@var{function_name}, @var{x}, @var{y}, @var{z})\n\
+Evaluate a function at a specific point of the domain and return the value. \n\
+With only two arguments, @var{Coordinate} is assumed to be an array holding \n\
+the coordinates of the point where the function should be evaluated.  With \n\
+three or more arguments, the coordinates are passed separately in @var{x}, \n\
+@var{y} and optionally @var{z}.  In the latter case and if the function \n\
+returns scalar values, @var{x}, @var{y} and @var{z} may also be arrays to \n\
+evaluate the function at multiple points at once.\n\
+@seealso{Function}\n\
+@end deftypefn")
+{
+
+  int nargin = args.length ();
+  octave_value retval = 0;
+  
+  if (nargin < 2)
+    print_usage ();
+  else
+    {
+      if (! function_type_loaded)
+        {
+          function::register_type ();
+          function_type_loaded = true;
+          mlock ();
+        }
+
+      if (args(0).type_id () == function::static_type_id ())
+        {
+          const function & fspo =
+            static_cast<const function&> (args(0).get_rep ());
+          const boost::shared_ptr<const dolfin::Function> 
+            & f = fspo.get_pfun ();
+          dim_vector dims;
+          dims.resize (2);
+          dims(0) = 1;
+          dims(1) = f->value_dimension (0);
+          Array<double> res(dims);
+
+          if (nargin == 2)
+            {
+              Array<double> point = args(1).array_value ();
+              if (!error_state)
+                {
+                  evaluate (*f, point, res);
+                  retval = octave_value (res);
+                }
+            }
+          else
+            {
+              dim_vector inDims;
+              inDims.resize (2);
+              inDims(0) = 1;
+              inDims(1) = nargin - 1;
+              Array<double> point(inDims);
+
+              std::vector<Array<double> > coords;
+              coords.reserve (inDims(1));
+              dim_vector argDims;
+              for (unsigned i = 1; i < nargin; ++i)
+                {
+                  coords.push_back (args(i).array_value ());
+                  const dim_vector curDims = coords.back ().dims ();
+                  if (i > 1 && argDims != curDims)
+                    {
+                      error ("feval: coordinate dimensions mismatch");
+                      break;
+                    }
+                  argDims = curDims;
+                }
+
+              if (res.nelem () != 1)
+                error ("feval: separate coordinate version only supported"
+                       "for scalar functions");
+
+              if (!error_state)
+                {
+                  Array<double> retValArray(argDims);
+
+                  for (size_t i = 0; i < retValArray.nelem (); ++i)
+                    {
+                      for (unsigned j = 0; j < coords.size (); ++j)
+                        point(j) = coords[j].fortran_vec ()[i];
+
+                      evaluate (*f, point, res);
+                      retValArray.fortran_vec ()[i] = res(0);
+                    }
+
+                  retval = octave_value (retValArray);
+                }
+            }
+        }
+    }
+
+  return retval;
+}