new file mode 100644
--- /dev/null
+++ b/include/interp_values.hpp
@@ -0,0 +1,144 @@
+#ifndef __INTERP_VALUES_HPP__
+#define __INTERP_VALUES_HPP__
+
+#include "linalg.hpp"
+#include <map>
+
+namespace bvp{
+  using std::map;
+  using namespace linalg;
+
+  class bdry_values;
+
+  //FIXME: Give this a better name
+  /// A unit normals class associated to an interpolator
+  class normals{
+  public:
+    /** \brief Return the kth coordinate of all normals along the
+    * boundary as interp values
+    */
+    bdry_values operator()(size_t k) const;
+    
+    /*! \name comparisons 
+     *
+     * Normals are the same if they correspond to the same domain
+     *(i.e. if their hash coincides)
+     */
+    //@{
+    bool operator!=(const normals& M) const;
+    bool operator==(const normals& M) const;
+    //@}
+  private:
+    template<typename RBF> friend class interpolator;
+    ///Only constructor available for interpolator class
+    normals(size_t rbfs_hash_in, const map<point, linalg::vector>&
+	    normals_in, size_t n_in);
+    ///Does nothing
+    normals(){};
+    ///Exception builder
+    badArgument not_initted() const;
+    
+    ///A hash used to return bdry_values objects
+    size_t rbfs_hash;
+
+    ///Number of interior points in the matrix, which interp_values need
+    size_t n;
+    ///Number of boundary points
+    size_t m;
+
+    ///Store normals in a matrix
+    linalg::matrix N;
+  };
+
+
+
+
+  /// 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;
+
+    interp_values operator+(const double a) const;
+    interp_values operator-(const double a) const;
+    interp_values operator*(const double a) const;
+    interp_values operator/(const double a) const;
+    //@}
+
+    /*! \name Non-member arithmetic operators
+     *
+     * For comfortable syntax.
+     */
+    //@{
+    friend interp_values operator+(double a, const interp_values& v); 
+    friend interp_values operator-(double a, const interp_values& v); 
+    friend interp_values operator*(double a, const interp_values& v); 
+    friend interp_values operator/(double a, const interp_values& v); 
+    //@}
+
+  protected:
+    /// Interpolator creates interp_values with this
+    interp_values(const linalg::vector& v_in, size_t rbfs_hash_in,
+		  size_t n_in, size_t m_in) 
+      : v(v_in), rbfs_hash(rbfs_hash_in), n(n_in), m(m_in) {};
+
+    /// No construction!
+    interp_values();
+
+  protected:
+    ///Exception builder
+    badArgument different_rbfs(int line, string file) const;
+
+    ///Interpolated values go here
+    linalg::vector v;
+    ///Two interp_values can be added iff these two values matches.
+    size_t rbfs_hash;
+    ///Number of interior points
+    size_t n;
+    ///Number of boundary points
+    size_t m;
+
+    template<typename RBF> friend class interpolator;
+  };
+  
+
+
+  /// For referring to values purely on the boundary of an interpolator.
+  class bdry_values : public interp_values{
+  public:
+    ///Empty initialisation
+    bdry_values(){};
+    ///Initialise the boundary values with interp_values
+    bdry_values(const interp_values& v_in);
+  private:
+    friend class normals;
+    ///Same as for parent class
+    bdry_values(const linalg::vector& v_in, size_t rbfs_hash_in,
+		size_t n_in, size_t m_in) 
+      : interp_values(v_in, rbfs_hash_in, n_in, m_in) {};
+  };
+
+
+
+  /*! \name Arithmetic operators
+   * For comfortable syntax.
+   */
+  //@{
+  interp_values operator+(double a, const interp_values& v); 
+  interp_values operator-(double a, const interp_values& v); 
+  interp_values operator*(double a, const interp_values& v); 
+  interp_values operator/(double a, const interp_values& v); 
+  //@}
+}
+
+#endif //__INTERP_VALUES_HPP__

new file mode 100644
--- /dev/null
+++ b/interp_values.cpp
@@ -0,0 +1,151 @@
+#include "include/interp_values.hpp"
+
+namespace bvp{
+
+  // ************************ interp_values stuff *****************
+  badArgument interp_values::different_rbfs(int line, string file) const
+  {
+    badArgument exc;
+    exc.reason = "Can't operate on interpolated values from "
+      "incompatible interpolators";
+    exc.line = line;
+    exc.file = file;
+    return exc;
+  }
+
+  interp_values 
+  interp_values::operator+(const interp_values& w) const
+  {
+    if(rbfs_hash != w.rbfs_hash)
+      throw different_rbfs(__LINE__, __FILE__);
+    
+    return interp_values(v + w.v,rbfs_hash,n,m);
+  }
+
+  interp_values 
+  interp_values::operator-(const interp_values& w) const
+  {
+    if(rbfs_hash != w.rbfs_hash)
+      throw different_rbfs(__LINE__, __FILE__);
+    
+    return interp_values(v - w.v,rbfs_hash,n,m);
+  }
+
+  interp_values 
+  interp_values::operator*(const interp_values& w) const
+  {
+    if(rbfs_hash != w.rbfs_hash)
+      throw different_rbfs(__LINE__, __FILE__);
+    return interp_values(v * w.v,rbfs_hash,n,m);
+  }
+
+  interp_values 
+  interp_values::operator/(const interp_values& w) const
+  {
+    if(rbfs_hash != w.rbfs_hash)
+      throw different_rbfs(__LINE__, __FILE__);
+    
+    return interp_values(v / w.v,rbfs_hash,n,m);
+  }
+
+  interp_values
+  interp_values::operator+(const double a) const
+  {
+    vector w(v.size(),a);
+    return interp_values(v+w,rbfs_hash,n,m);
+  }
+  
+  interp_values
+  interp_values::operator-(const double a) const
+  {
+    vector w(v.size(),a);
+    return interp_values(v-w, rbfs_hash,n,m);
+  }
+
+  interp_values
+  interp_values::operator*(const double a) const
+  {
+    return interp_values(v*a,rbfs_hash,n,m);
+  }
+  
+  interp_values
+  interp_values::operator/(const double a) const
+  {
+    return interp_values(v/a, rbfs_hash,n,m);
+  }
+
+  //*************** bdry_vals stuff ************************************
+
+  bdry_values::bdry_values(const interp_values& v_in) 
+  : interp_values(v_in)
+  {
+    if(v.size() == m+n){       
+      //Lower part of the vector contains the boundary
+      linalg::slice s(n+1,v.size());
+      v = v(s);
+    }
+    else if(v.size() != m){ 
+      throw badArgument("Cannot initialise bdry_values from given "
+			"interp_values", __FILE__,__LINE__);
+    }
+  }
+
+  //*************** normals stuff *************************************
+  
+  normals::normals(size_t rbfs_hash_in, 
+	       const map<point, linalg::vector>& normals_in, 
+	       size_t n_in) : rbfs_hash(rbfs_hash_in), n(n_in)
+  {
+    map<point, linalg::vector>::const_iterator I;
+    size_t rows = normals_in.size();
+    size_t cols = normals_in.begin() -> second.size();
+    slice s(1,cols);
+    size_t i = 1;
+    matrix N_temp(rows,cols);
+    for(I = normals_in.begin(); I != normals_in.end(); I++){
+      N_temp(i,s) = (I->second)/(norm(I->second));
+      i++;
+    }
+    N = N_temp;
+  }
+
+  bdry_values normals::operator()(size_t k) const
+  {
+    slice s(1,N.rows());
+    return bdry_values(N(s,k),rbfs_hash,n,m);
+  }
+
+  bool normals::operator==(const normals& M) const
+  {
+    return rbfs_hash == M.rbfs_hash;
+  }
+
+  bool normals::operator!=(const normals& M) const
+  {
+    return rbfs_hash != M.rbfs_hash;
+  }
+
+  //*************** Non-member global friend functions ******************
+
+  interp_values operator+(double a, const interp_values&  v)
+  { 
+    return v+a;
+  } 
+
+  interp_values operator-(double a, const interp_values&  v)
+  { 
+    linalg::vector w(v.v.size(),a);
+    return interp_values(w-v.v, v.rbfs_hash,v.n,v.m);
+  } 
+
+  interp_values operator*(double a, const interp_values&  v)
+  { 
+    return v*a;
+  } 
+
+  interp_values operator/(double a, const interp_values&  v)
+  { 
+    linalg::vector w(v.v.size(),a);
+    return interp_values(w/v.v, v.rbfs_hash,v.n,v.m);
+  } 
+}