--- a/include/interpolator.hpp
+++ b/include/interpolator.hpp
@@ -13,6 +13,8 @@
   using std::map;
   using boost::shared_ptr;
 
+  class interp_values;
+
   template<typename RBF>
   class interpolator : public realfunc{
   public:
@@ -49,20 +51,38 @@
     //@{
     void interpolate(const map<point, double>& Xi); 
     void interpolate(shared_ptr<linear_BVP2> bvp);
+    
+    /** Interpolate given values on the domain
+     * 
+     * The given values must match the values returned by other
+     * interpolators on compatible domains.
+     */
+    void interpolate(const interp_values& values); 
     //@}
     
     /** @name Evaluations and derivatives
+     * 
+     * Full domain evaluations and derivatives are 
      */
     //@{
-    ///Evaluation
+    ///Point evaluation
     double operator()(const point& p) const;
-    ///Evaluation
+    ///Point evaluation
     double at(const point& p) const;
 
-    /// First derivative
+    ///Full domain evaluation
+    interp_values operator()() const; 
+    ///Full domain evaluation
+    interp_values at() const; 
+
+    /// Pointwise first derivative
     double d(const point& p, size_t k) const;
-    /// Second derivatives
+    /// Full domain first erivative
+    interp_values d(size_t k) const; 
+    /// Pointwise second derivatives
     double d2(const point &p, size_t k1, size_t k2) const;
+    /// Full domain second derivatives
+    interp_values d2(size_t k1, size_t k2) const; 
     //@}
 
     /** @name Precomputation of RBFs
@@ -130,7 +150,13 @@
 
     ///Is the interpolator ready for use?
     bool initted;
-    void not_initted(int line, string file) const; //Exception thrower
+
+    /** @name Exception throwers
+     */
+    //@{
+    void not_initted(int line, string file) const; 
+    void not_precomputed(int line, string file) const;
+    //@}
     
     ///Coefficients of the RBFs
     linalg::vector coeffs;
@@ -150,26 +176,67 @@
      * evaluations of the interpolator when the interpolator's domain
      * doesn't change.
      *
-     * The keys in this vector are vectors that represent the
-     * multindex for the derivative, where missing trailing entries in
-     * the vector represents zeros. Thus, an empty vector represents
+     * The keys in this map are vectors that represent the multindex
+     * for the derivative, where missing trailing entries in the
+     * vector represents zeros. Thus, an empty vector represents
      * evaluation instead of derivatives.
      */
     mutable map<std::vector<size_t>, matrix> precomp_rbfs;
     
-    /// Have the RBFs or their derivatives been precomputed?
+    /** @name  Evaluation flags
+     *Have the RBFs or their derivatives been precomputed?
+     */
+    //@{
     bool ev_precomp;
     bool d1_precomp;
     bool d2_precomp;
+    //@}
 
     /// Precomputed values using precomp_rbfs
     mutable map<std::vector<size_t>, map<point, double> >
     precomp_values;
+    /// Same, but storing the vectors.
+    mutable map<std::vector<size_t>, linalg::vector >
+    precomp_values_vec;
     /// Invert phis\coeffs and return result vector in a map.
     map<point, double> precompute_values(const matrix& phis) const;
   };
 
-  //For comfortable syntax
+  /// Interpolated values manipulated by the interpolator class.
+  class interp_values{
+  public:
+    /** @name Arithmetic with interpolated values
+     *
+     * These operators allow pointwise arithmetic with interpolated
+     * values returned by interpolators. They must work with
+     * interpolated values from compatible interpolators
+     * (i.e. interpolators with the same RBFs).
+     */
+    //@{
+    interp_values operator+(const interp_values& w) const;
+    interp_values operator-(const interp_values& w) const;
+    interp_values operator*(const interp_values& w) const;
+    interp_values operator/(const interp_values& w) const;
+    //@}
+  private:
+    
+    /// No construction!
+    interp_values();
+    /// Interpolator creates interp_values with this
+    interp_values(const linalg::vector& v_in, size_t rbfs_hash_in) 
+      : v(v_in), rbfs_hash(rbfs_hash_in) {};
+
+    ///Interpolated values go here
+    linalg::vector v;
+    ///Two interp_values can be added iff these two values matches.
+    size_t rbfs_hash;
+    ///Exception thrower
+    void different_rbfs(int line, string file) const;
+
+    template<typename RBF> friend class interpolator;
+  };
+
+  ///For comfortable syntax
   template <typename RBF>
   interpolator<RBF> operator*(double a, const interpolator<RBF>& u)
   {