new file mode 100644
--- /dev/null
+++ b/.hgignore
@@ -0,0 +1,7 @@
+syntax: glob
+html/*
+*.ogg
+*.ogv
+*.matrix
+*.map
+*.o

--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,6 @@
 CPP = g++
 LINKING = -lgsl -lgslcblas 
-CFLAGS = -ggdb3
+CFLAGS = -O3
 OPTIONS = -Wall -pedantic  -W -Werror -Wconversion -Wshadow \
 	  -Wpointer-arith -Wcast-qual -Wcast-align -Wwrite-strings  \
 	  -fshort-enums -fno-common -Wfatal-errors 

new file mode 100644
--- /dev/null
+++ b/data/linear_wave.m
@@ -0,0 +1,64 @@
+n = 20
+m = 6;
+r = linspace(0,1,n);
+
+rrintr = 0;
+ththintr = 0;
+for k = 2:n-1
+  c = (k-1)*m;
+  th = [0:2*pi/c:2*pi*(1-1/c)];
+  for j = 1:c
+    rrintr(end+1) = r(k);
+    ththintr(end+1) = th(j);
+  endfor
+endfor
+#Remove one in order to get a nice non-prime number.
+#rrintr = rrintr(2:end);
+#ththintr = ththintr(2:end);
+
+r=1;
+c = (n-1)*m;
+th = [0:2*pi/c:2*pi*(1-1/c)];
+rrbdry = 1;
+ththbdry = 0;
+for j = 2:c
+  rrbdry(end+1) = r;
+  ththbdry(end+1) = th(j);
+endfor
+[xxintr,yyintr] = pol2cart(ththintr, rrintr);
+[xxbdry,yybdry] = pol2cart(ththbdry, rrbdry);
+intr = [xxintr',yyintr'];
+bdry = [xxbdry',yybdry'];
+
+##gotta avoid zero coordinates, so shift the boundary slightly
+bdry = ([cos(pi/200), sin(pi/200); -sin(pi/200) cos(pi/200)]*bdry')';
+nrml = [bdry, bdry];
+save "circ_intr.matrix" intr;
+save "circ_bdry.matrix" bdry;
+save "circ_nrml.matrix" nrml;
+
+xxbdry = bdry(:,1)';
+yybdry = bdry(:,2)';
+
+##radial
+##u0 = (1-(xxintr.^2 + yyintr.^2)')/2; 
+
+##nonradial
+h0 = [];
+ for i = 1:columns(xxintr)
+   h0(i) = fr( (xxintr(i)-0.25)^2 + yyintr(i)^2, 0.025);
+ endfor
+ h0 = 4*h0(:)+0.1;
+
+h0 = [h0];
+
+h_init = [xxintr', yyintr', h0];
+
+save "wave_init.map" h_init;
+
+
+hold off;
+plot(intr(:,1),intr(:,2),'xr');
+ hold on;
+plot(bdry(:,1),bdry(:,2),'xb');
+size([intr;bdry])

--- a/include/interp_values.hpp
+++ b/include/interp_values.hpp
@@ -15,8 +15,8 @@
   class normals{
   public:
     /** \brief Return the kth coordinate of all normals along the
-    * boundary as interp values
-    */
+     * boundary as interp values
+     */
     bdry_values operator()(size_t k) const;
     
     /*! \name comparisons 
@@ -32,7 +32,7 @@
     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);
+            normals_in, size_t n_in);
     ///Does nothing
     normals(){};
     ///Exception builder
@@ -89,7 +89,7 @@
   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) 
+                  size_t n_in, size_t m_in) 
       : v(v_in), rbfs_hash(rbfs_hash_in), n(n_in), m(m_in) {};
 
     /// No construction!
@@ -124,10 +124,26 @@
     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) 
+                size_t n_in, size_t m_in) 
       : interp_values(v_in, rbfs_hash_in, n_in, m_in) {};
   };
 
+  /// For referring to values purely on the interior of an interpolator.
+  class intr_values : public interp_values{
+  public:
+    ///Empty initialisation
+    intr_values(){};
+    ///Initialise the interior values with interp_values
+    intr_values(const interp_values& v_in);
+  private:
+    ///Same as for parent class
+    intr_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

--- a/include/interpolator.hpp
+++ b/include/interpolator.hpp
@@ -95,7 +95,6 @@
     void interpolate(const interp_values& values); 
     //@}
     
-
     /** @name Pointwise evaluation and derivatives
      */
     //@{
@@ -122,7 +121,6 @@
     interp_values d2(size_t k1, size_t k2) const; 
     //@}
 
-
     /** @name Precomputation of RBFs
      *
      * The following functions will precompute upon request all the
@@ -146,6 +144,8 @@
      * and for other methods. They do not factor a matrix again.
      */
     //@{
+    void set_f(const intr_values &f);
+    void set_g(const bdry_values &g);
     void set_f(const realfunc &f);
     void set_g(const realfunc &g);
     void set_f(const map<point, double>& f);
@@ -183,9 +183,18 @@
      */
     normals get_normals() const;
     
+    /** \brief This will put into an ostream a linalg::matrix where
+     *   the first \f$n-1\f$ columns will be all the points where the
+     *   interpolator is defined and the \f$n^{th}\f$ column is the
+     *   values that the interpolator takes on that point.
+     */
+    friend
+    template<typename RBF>
+    linalg::matrix& operator<< (ostream& os, 
+    
   private:
     ///Once the matrix is defined, this function inverts it.
-    void computecoeffs(); 
+    void computecoeffs(bool rhs_defined = false); 
 	 
     ///Perform the actual interpolation.
     void init(shared_ptr<linear_BVP2> bvp);
@@ -207,6 +216,7 @@
     ///Is the interpolator ready for use?
     bool initted;
 
+
     /** @name Exception builders
      */
     //@{

--- a/include/linalg.hpp
+++ b/include/linalg.hpp
@@ -18,6 +18,8 @@
 #include <gsl/gsl_permutation.h>
 #include "error.hpp"
 
+#include <boost/operators.hpp>
+
 /// Linear algebra namespace
 namespace linalg{
   using namespace error_handling;

--- a/interp_values.cpp
+++ b/interp_values.cpp
@@ -79,7 +79,7 @@
   //*************** bdry_vals stuff ************************************
 
   bdry_values::bdry_values(const interp_values& v_in) 
-  : interp_values(v_in)
+    : interp_values(v_in)
   {
     if(v.size() == m+n){       
       //Lower part of the vector contains the boundary
@@ -88,23 +88,37 @@
     }
     else if(v.size() != m){ 
       throw badArgument("Cannot initialise bdry_values from given "
-			"interp_values", __FILE__,__LINE__);
+                        "interp_values", __FILE__,__LINE__);
     }
     // Else, we're initting with interp_values that are the right size
     // and there's no need to do anything, since they're bona fide
     // boundary values (parent copy ctor does the work).
+  }
 
-      //debug
-      using namespace std; 
-      cout << "The boundary vector: " << endl
-	   << v << endl;
+  //*************** intr_vals stuff *****************************
+
+  intr_values::intr_values(const interp_values& v_in) 
+    : interp_values(v_in)
+  {
+    if(v.size() == m+n){       
+      //Upper part of the vector contains the boundary
+      linalg::slice s(1,n);
+      v = v(s);
+    }
+    else if(v.size() != n){ 
+      throw badArgument("Cannot initialise intr_values from given "
+                        "interp_values", __FILE__,__LINE__);
+    }
+    // Else, we're initting with interp_values that are the right size
+    // and there's no need to do anything, since they're bona fide
+    // boundary values (parent copy ctor does the work).
   }
 
   //*************** 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)
+                   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();

--- a/interpolator.cpp
+++ b/interpolator.cpp
@@ -33,20 +33,20 @@
 
   template<typename RBF>
   interpolator<RBF>::interpolator(shared_ptr<domain> Omega, 
-				  const map<point, double>& Xi)
+                                  const map<point, double>& Xi)
   {
     //Check that Xi matches the domain in size
     size_t Omega_size = Omega -> get_interior().size() 
-                       + Omega -> get_boundary().size();
+      + Omega -> get_boundary().size();
     if(Xi.size() != Omega_size)
     {
       badArgument exc;
       if(Xi.size() < Omega_size)
-	exc.reason = "Did not provide enough interpolation data for every " 
-	  "point in the given domain.";
+        exc.reason = "Did not provide enough interpolation data for every " 
+          "point in the given domain.";
       else
-	exc.reason = "Provided more interpolation data than points in the " 
-	  "given domain.";
+        exc.reason = "Provided more interpolation data than points in the " 
+          "given domain.";
       exc.line = __LINE__;
       exc.file = __FILE__;
       throw exc;
@@ -59,19 +59,19 @@
 
     //Extract f and g from given information
     for(map<point, double>::const_iterator I = Xi.begin();
-	I != Xi.end(); I++)
+        I != Xi.end(); I++)
     {
       if(utils::contains(Omega -> get_interior(), I -> first))
-	f[I -> first] = I -> second;
+        f[I -> first] = I -> second;
       else if(utils::contains(Omega -> get_boundary(), I -> first))
-	g[I -> first] = I -> second;
+        g[I -> first] = I -> second;
       else
       {
-	badArgument exc;
-	exc.reason = "The interpolation data contains points not in the given domain.";
-	exc.line = __LINE__;
-	exc.file = __FILE__;
-	throw exc;
+        badArgument exc;
+        exc.reason = "The interpolation data contains points not in the given domain.";
+        exc.line = __LINE__;
+        exc.file = __FILE__;
+        throw exc;
       }
     } 
     
@@ -110,15 +110,15 @@
     size_t dimension = 0;
     for(map<point,double>::const_iterator i = Xi.begin(); i != Xi.end(); i++){
       if(!dim_set){
-	dimension = (i -> first).size();
-	dim_set = true;
+        dimension = (i -> first).size();
+        dim_set = true;
       }
       else if(dimension != (i -> first).size()){
-	badArgument exc;
-	exc.reason = "Inconformant dimensions in interpolation data.";
-	exc.line = __LINE__;
-	exc.file = __FILE__;
-	throw exc;
+        badArgument exc;
+        exc.reason = "Inconformant dimensions in interpolation data.";
+        exc.line = __LINE__;
+        exc.file = __FILE__;
+        throw exc;
       }
       intr.insert( i->first);
     }
@@ -140,7 +140,7 @@
     if(!initted){
       badArgument exc;
       exc.reason = "Cannot interpolate with other interpolator's "
-	"values without defining the domain.";
+        "values without defining the domain.";
       exc.line = __LINE__;
       exc.file = __FILE__;
       throw exc;
@@ -148,7 +148,7 @@
     if(rbfs_hash != values.rbfs_hash){
       badArgument exc;
       exc.reason = "Cannot interpolate with values from incompatible "
-	"interpolator." ;
+        "interpolator." ;
       exc.line = __LINE__;
       exc.file = __FILE__;
       throw exc;
@@ -202,7 +202,7 @@
     I = interior.begin();
     for(size_t i = 1; i <= n; i++){
       for(size_t j = 1; j <= n+m; j++)
-	M(i,j) = L -> at(rbfs[j-1], *I);
+        M(i,j) = L -> at(rbfs[j-1], *I);
       I++;
     }
     
@@ -210,7 +210,7 @@
     J = normals.begin();
     for(size_t i = n+1; i <= n+m; i++){
       for(size_t j = 1; j <= n+m; j++)
-	M(i,j) = B -> at(rbfs[j-1], J->first, J->second);
+        M(i,j) = B -> at(rbfs[j-1], J->first, J->second);
       J++;
     }
    
@@ -258,14 +258,14 @@
       set<point>::iterator I;
       size_t i;
       for(I = thebvp -> get_domain () -> get_interior().begin(), i = 1; 
-	  i <= n; i++, I++)
-	for(size_t j = 0; j < n+m; j++)
-	  phis(i,j+1) = rbfs[j].at(*I);
+          i <= n; i++, I++)
+        for(size_t j = 0; j < n+m; j++)
+          phis(i,j+1) = rbfs[j].at(*I);
       
       for(I = thebvp -> get_domain() -> get_boundary().begin(), i=n+1; 
-	  i <= n+m; i++, I++)
-	for(size_t j = 0; j < n+m; j++)
-	  phis(i,j+1) = rbfs[j].at(*I);
+          i <= n+m; i++, I++)
+        for(size_t j = 0; j < n+m; j++)
+          phis(i,j+1) = rbfs[j].at(*I);
       
       precomp_rbfs[alpha] = phis;
       computecoeffs();
@@ -286,25 +286,25 @@
       //No, so precompute the RBFs first derivatives at all points of
       //the domain.
       for(size_t k = 1; k <= thebvp -> get_domain() -> get_dimension(); 
-	  k++){
-	std::vector<size_t> alpha(k); alpha[k-1]++;
-	matrix phis(n+m,n+m);
-	{
-	  set<point>::iterator I;
-	  size_t i;
-	  for(I = thebvp -> get_domain() -> get_interior().begin(), 
-		i = 1; i <= n; i++, I++)
-	    for(size_t j = 0; j < n+m; j++)
-	      phis(i,j+1) = rbfs[j].d(*I,k);
+          k++){
+        std::vector<size_t> alpha(k); alpha[k-1]++;
+        matrix phis(n+m,n+m);
+        {
+          set<point>::iterator I;
+          size_t i;
+          for(I = thebvp -> get_domain() -> get_interior().begin(), 
+                i = 1; i <= n; i++, I++)
+            for(size_t j = 0; j < n+m; j++)
+              phis(i,j+1) = rbfs[j].d(*I,k);
       
-	  for(I = thebvp -> get_domain() -> get_boundary().begin(), 
-		i=n+1; i <= n+m; i++, I++)
-	    for(size_t j = 0; j < n+m; j++)
-	      phis(i,j+1) = rbfs[j].d(*I,k);
-	}
+          for(I = thebvp -> get_domain() -> get_boundary().begin(), 
+                i=n+1; i <= n+m; i++, I++)
+            for(size_t j = 0; j < n+m; j++)
+              phis(i,j+1) = rbfs[j].d(*I,k);
+        }
 
-	precomp_rbfs[alpha] = phis;
-	computecoeffs();
+        precomp_rbfs[alpha] = phis;
+        computecoeffs();
       }
     }
     d1_precomp = true;
@@ -319,24 +319,24 @@
       size_t dim = thebvp -> get_domain() -> get_dimension();
       //precompute all second derivatives
       for(size_t k1 = 1; k1 <= dim; k1++){
-	for(size_t k2 = k1; k2 <= dim; k2++){
-	  std::vector<size_t> alpha(k1>k2?k1:k2); alpha[k1-1]++; alpha[k2-1]++;
-	    matrix phis(n+m,n+m);
-	    set<point>::iterator I;
-	    size_t i;
-	    for(I = thebvp -> get_domain () -> get_interior().begin(), i = 1; 
-		i <= n; i++, I++)
-	      for(size_t j = 0; j < n+m; j++)
-		phis(i,j+1) = rbfs[j].d2(*I,k1,k2);
+        for(size_t k2 = k1; k2 <= dim; k2++){
+          std::vector<size_t> alpha(k1>k2?k1:k2); alpha[k1-1]++; alpha[k2-1]++;
+          matrix phis(n+m,n+m);
+          set<point>::iterator I;
+          size_t i;
+          for(I = thebvp -> get_domain () -> get_interior().begin(), i = 1; 
+              i <= n; i++, I++)
+            for(size_t j = 0; j < n+m; j++)
+              phis(i,j+1) = rbfs[j].d2(*I,k1,k2);
       
-	    for(I = thebvp -> get_domain() -> get_boundary().begin(), i=n+1; 
-		i <= n+m; i++, I++)
-	      for(size_t j = 0; j < n+m; j++)
-		phis(i,j+1) = rbfs[j].d2(*I,k1,k2);
+          for(I = thebvp -> get_domain() -> get_boundary().begin(), i=n+1; 
+              i <= n+m; i++, I++)
+            for(size_t j = 0; j < n+m; j++)
+              phis(i,j+1) = rbfs[j].d2(*I,k1,k2);
       
-	    precomp_rbfs[alpha] = phis;
-	    computecoeffs();	  
-	}
+          precomp_rbfs[alpha] = phis;
+          computecoeffs();	  
+        }
       }
     }
     d2_precomp = true;
@@ -470,7 +470,7 @@
     if(this -> rbfs_hash != u.rbfs_hash){
       badArgument exc;
       exc.reason = 
-    	"Cannot add interpolators on different domains (interior and boundary must match)."; 
+        "Cannot add interpolators on different domains (interior and boundary must match)."; 
       exc.line = __LINE__;
       exc.file = __FILE__;
       throw exc;	
@@ -490,7 +490,7 @@
     if(this -> rbfs_hash != u.rbfs_hash){
       badArgument exc;
       exc.reason = 
-    	"Cannot subtract interpolators on different domains (interior and boundary must match)."; 
+        "Cannot subtract interpolators on different domains (interior and boundary must match)."; 
       exc.line = __LINE__;
       exc.file = __FILE__;
       throw exc;	
@@ -525,7 +525,7 @@
   normals interpolator<RBF>::get_normals() const{
     if(!initted){
       throw badArgument("Interpolator not initialised, so there are "
-			"no normals!",__FILE__,__LINE__);
+                        "no normals!",__FILE__,__LINE__);
     }
     return nrmls;
   }
@@ -533,6 +533,34 @@
   //******************* data manipulation **************************
 
   template<typename RBF>
+  void interpolator<RBF>::set_f(const intr_values& f){
+    if(!initted){
+      throw not_initted(__LINE__, __FILE__);
+    }
+
+    std::vector<size_t> alpha;
+    slice s(1,n);
+    precomp_values_vec[alpha](s) = f.v;
+    //precomp_values[alpha] = valsvec2map(precomp_values_vec[alpha]);
+
+    computecoeffs(true);
+  }
+
+  template<typename RBF>
+  void interpolator<RBF>::set_g(const bdry_values& g){
+    if(!initted){
+      throw not_initted(__LINE__, __FILE__);
+    }
+    
+    std::vector<size_t> alpha;
+    slice s(n+1,n+m);
+    precomp_values_vec[alpha](s) = g.v;
+    //precomp_values[alpha] = valsvec2map(precomp_values_vec[alpha]);
+    
+    computecoeffs(true);
+  }
+
+  template<typename RBF>
   void interpolator<RBF>::set_f(const realfunc& f){
     if(!initted){
       throw not_initted(__LINE__, __FILE__);
@@ -611,23 +639,28 @@
   }
 
   template<typename RBF>
-  void interpolator<RBF>::computecoeffs()
+  void interpolator<RBF>::computecoeffs(bool rhs_defined)
   {
     linalg::vector rhs(n+m);
 
     //Compute the RBF coefficients
+    if(rhs_defined){
+      std::vector<size_t> alpha;
+      rhs = precomp_values_vec[alpha];
+    }
+    else
     {
       map<point, double>::const_iterator I;
 
       I = (thebvp -> get_f()).begin();
       for(size_t i = 1; i <= n; i++){
-	rhs(i) = I->second;
-	I++;
+        rhs(i) = I->second;
+        I++;
       }
       I = (thebvp -> get_g()).begin();   
       for(size_t i = n+1; i <= n+m; i++){
-	rhs(i) = I->second;
-	I++;
+        rhs(i) = I->second;
+        I++;
       }
     }
     coeffs = M.inv(rhs);
@@ -638,9 +671,9 @@
       typename map<std::vector<size_t>, matrix>::const_iterator I;
       
       for(I = precomp_rbfs.begin(); I != precomp_rbfs.end(); I++){
-	linalg::vector vals = (I->second)*coeffs;
-	precomp_values_vec[I -> first] = vals;
-	precomp_values[I -> first] = valsvec2map(vals);
+        linalg::vector vals = (I->second)*coeffs;
+        precomp_values_vec[I -> first] = vals;
+        precomp_values[I -> first] = valsvec2map(vals);
       }
     }
   }
@@ -654,12 +687,12 @@
       set<point>::const_iterator I;
       size_t i;
       for(I = thebvp -> get_domain() -> get_interior().begin(), i = 1;
-	  i <= n; I++, i++)
-	out[*I] = values(i);
+          i <= n; I++, i++)
+        out[*I] = values(i);
       
       for(I = thebvp -> get_domain() -> get_boundary().begin(), i = n+1;
-	  i <= n+m; I++, i++)
-	out[*I] = values(i);
+          i <= n+m; I++, i++)
+        out[*I] = values(i);
     }
     return out;
   }
@@ -680,11 +713,11 @@
 
 /*! \mainpage Outline of the RBF-DDM code
 
-Thank you for your interest in my code. Above you can browse the
-structure of my code. You can also read here a technical design
-document that also outlines certain difficulties I have been facing
-with it:
+  Thank you for your interest in my code. Above you can browse the
+  structure of my code. You can also read here a technical design
+  document that also outlines certain difficulties I have been facing
+  with it:
 
-http://inversethought.com/jordi/design.pdf
+  http://inversethought.com/jordi/design.pdf
 */
 

--- a/main-sw-rk4.cpp
+++ b/main-sw-rk4.cpp
@@ -1,4 +1,5 @@
 #include <iostream>
+#include <iomanip>
 #include <fstream>
 #include <map>
 #include <set>
@@ -145,7 +146,10 @@
     std::vector<interpolator<RBF_TYPE> >
       // FIXME: make the copy ctor for interps make different copies
       // of the bvp.
-      k1(3,h0), k2(3,h0), k3(3,h0), k4(3,h0); 
+      k1(3,u0), k2(3,u0), k3(3,u0), k4(3,u0); 
+
+    //debug
+    cout << "Size of interpolator is " << sizeof(k1[0]) << endl;
 
     Fu<RBF_TYPE> f_u;
     Fv<RBF_TYPE> f_v;
@@ -159,13 +163,13 @@
     f_h.set_dt(dt);
 
     //main loop, timestepping with RK4.
-    size_t maxiter = 5;
+    size_t maxiter = 1000;
     for(size_t i = 1; i <= maxiter; i++){ 
       cout << "Now on iteration #" << i << endl;
-      if(i % 1 == 0){
- 	save_interp(Omega,u0,i,"u");
- 	save_interp(Omega,v0,i,"v");
- 	save_interp(Omega,h0,i,"h");
+      if(i % 10 == 0){
+        save_interp(Omega,u0,i,"u");
+        save_interp(Omega,v0,i,"v");
+        save_interp(Omega,h0,i,"h");
       }
 
       //k1 
@@ -180,10 +184,7 @@
       set_bdry_conds(k1[0],k1[1]);
       
       cout << "k1 done!" << endl;
-
-      //debug Euler
-      u0 = k1[0]; v0 = k1[1]; h0 = k1[2]; continue;
-           
+          
       //k2
       f_u.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
       f_v.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
@@ -301,25 +302,14 @@
 		 const interpolator<RBF>& u, 
 		 size_t n, string which_interp)
 {
-
-
   using namespace std;
   string filename = "results/" + which_interp;
-  if(n < 10000)
-    filename += "0";
-  if(n < 1000)
-    filename += "0";
-  if(n < 100)
-    filename += "0";
-  if(n < 10)
-    filename += "0";
   stringstream ss;
+  ss << std::setw(5) << setfill('0') << n;
   string n_string;
-  ss << n;
   ss >> n_string;
   filename += n_string + ".map";
-  ofstream ofs(filename.c_str());
-    
+  ofstream ofs(filename.c_str());    
   {
     slice s(1,2);
     matrix M(Omega -> get_interior().size() +
@@ -375,7 +365,9 @@
   //Project the vector (u,v) onto the perp of (N(1),N(2)) = (-N(2),N(1)).
   bvp::bdry_values u_bdry_new = u_bdry*N(2)*N(2) - v_bdry*N(1)*N(2);
   bvp::bdry_values v_bdry_new = v_bdry*N(1)*N(1) - u_bdry*N(1)*N(2);
-  
+
   u.set_bdry_values(u_bdry_new);
   v.set_bdry_values(v_bdry_new);
 }
+
+

--- a/main.cpp
+++ b/main.cpp
@@ -21,87 +21,25 @@
 
 #define RBF_TYPE rbf::conical
 
-//Will solve the system
-//
-//    u_t = -u*u_x - v*u_y - g*h_x         =: f1(u,v,h)
-// 
-//    v_t = -u*v_x - v*v_y - g*h_y         =: f2(u,v,h)
-//
-//    h_t = -u*h_x - h*u_x - v*h_y - h*v_y =: f3(u,v,h)
-//
-//where g is the acceleration due to gravity, u(x,y,0) = v(x,y,0) = 0,
-//and h starts out as an undisturbed surface high with a small
-//disturbance. We'll do this by RK4.
-
-
-const double g = 9.8; // m/s^2
-
 const double pi = 3.141592653589793238462643383279502;
 
-template<typename RBF>
-class Fgen {
+class wave_op : public linear_diff_op2{
 public:
-  void set_interps(const interpolator<RBF>& u0, 
-		   const interpolator<RBF>& v0,
-		   const interpolator<RBF>& h0);
-  void set_u(const interpolator<RBF>& u0);
-  void set_v(const interpolator<RBF>& v0);
-  void set_h(const interpolator<RBF>& h0);
-  void set_dt(double dt_in);
-protected:
-  interpolator<RBF> u;
-  interpolator<RBF> v;
-  interpolator<RBF> h;
-  double dt;
-};
-
-template<typename RBF>
-class Fu : public Fgen<RBF>{
-  using Fgen<RBF>::u;
-  using Fgen<RBF>::v;
-  using Fgen<RBF>::h;
-  using Fgen<RBF>::dt;
-public:
-  bvp::interp_values at() const;
-};
-
-template<typename RBF>
-class Fv : public Fgen<RBF>{
-  using Fgen<RBF>::u;
-  using Fgen<RBF>::v;
-  using Fgen<RBF>::h;
-  using Fgen<RBF>::dt;
-public:
-  bvp::interp_values at() const;
-};
-
-template<typename RBF>
-class Fh : public Fgen<RBF>{
-  using Fgen<RBF>::u;
-  using Fgen<RBF>::v;
-  using Fgen<RBF>::h;
-  using Fgen<RBF>::dt;
-public:
-  bvp::interp_values at() const;
+  wave_op(double dt_in): dt2(dt_in*dt_in){};
+  double at(const realfunc& f, const point& p) const
+  {
+    return dt2*L(f,p) + f(p);
+  }
+private:
+  Laplacian L;
+  double dt2;
 };
 
 class zero_func : public realfunc{
 public:
-  double at(const point& p) const{p.size(); return 0;};
+  double at(const point&) const{return 0;};
 };
 
-//************ Function declarations ****************************
-
-template<typename RBF>
-void save_interp(const shared_ptr<domain> Omega, 
-		 const interpolator<RBF>& u, 
-		 size_t n, string which_interp);
-
-template<typename RBF>
-void set_bdry_conds(interpolator<RBF> &u, interpolator<RBF> &v);
-
-
-//********************** Main ******************************************
 int main()
 {  
   gsl_set_error_handler(&error_handling::errorHandler);
@@ -110,137 +48,37 @@
     using namespace std;
     using boost::shared_ptr;
    
-    map<point, double> h_init = utils::read_pd_map("data/h_init.map");
+    map<point, double> u_init = utils::read_pd_map("data/wave_init.map");
       
     shared_ptr<domain> Omega(new domain("data/circ_intr.matrix",
-					"data/circ_bdry.matrix",
-					"data/circ_nrml.matrix"));
-
-    shared_ptr<Id_op> I(new Id_op);
-    shared_ptr<dirichlet_op> D(new dirichlet_op);
-    zero_func Z;
-    //FIXME: Turns this into a purely interpolating interpolator,
-    //without a BVP.
-    shared_ptr<linear_BVP2> 
-      u_bvp_init(new linear_BVP2(Omega, I, D, Z, Z)),
-      v_bvp_init(new linear_BVP2(Omega, I, D, Z, Z));
-      
+                                        "data/circ_bdry.matrix",
+                                        "data/circ_nrml.matrix"));
     
-    //init the interps.
-    using namespace rbf;
-    conical::set_n(7); thin_plate_spline::set_n(6); 
-    c_infty_rbf::set_epsilon(0.01);
-
-    cout << "Initialising... please wait." << endl;
-
-    interpolator<RBF_TYPE> 
-      u0(u_bvp_init), 
-      v0(v_bvp_init), 
-      h0(Omega, h_init);
-
-    //Precompute some RBFs
-    u0.precompute_ev(); v0.precompute_ev(), h0.precompute_ev();
-    u0.precompute_d1(); v0.precompute_d1(), h0.precompute_d1();
+    zero_func g;
 
-    //Intermediate interpolators for RK4
-    std::vector<interpolator<RBF_TYPE> >
-      // FIXME: make the copy ctor for interps make different copies
-      // of the bvp.
-      k1(3,u0), k2(3,u0), k3(3,u0), k4(3,u0); 
-
-    //debug
-    cout << "Size of interpolator is " << sizeof(k1[0]) << endl;
-
-    Fu<RBF_TYPE> f_u;
-    Fv<RBF_TYPE> f_v;
-    Fh<RBF_TYPE> f_h;
-
-    //Negative because of the way the F's below are written.
-    double dt = -0.01;
-
-    f_u.set_dt(dt);
-    f_v.set_dt(dt);
-    f_h.set_dt(dt);
-
-    //main loop, timestepping with RK4.
-    size_t maxiter = 1000;
-    for(size_t i = 1; i <= maxiter; i++){ 
-      cout << "Now on iteration #" << i << endl;
-      if(i % 1 == 0){
- 	save_interp(Omega,u0,i,"u");
- 	save_interp(Omega,v0,i,"v");
- 	save_interp(Omega,h0,i,"h");
-      }
+    //timestep
+    double dt=1e-2;
+    
+    shared_ptr<wave_op> W(new wave_op(dt) );
+    shared_ptr<dirichlet_op> D(new dirichlet_op);
+    shared_ptr<linear_BVP2> the_bvp(new linear_BVP2(Omega, W, D,
+                                                    u_init, g));
 
-      //k1 
-      f_u.set_interps(u0,v0,h0);
-      f_v.set_interps(u0,v0,h0);
-      f_h.set_interps(u0,v0,h0);
-
-//       k1[0].interpolate(f_u.at());
-//       k1[1].interpolate(f_v.at());
-//       k1[2].interpolate(f_h.at());
-
-      //debug Euler
-      k1[0].interpolate(u0.at() + f_u.at());
-      k1[1].interpolate(v0.at() + f_v.at());
-      k1[2].interpolate(h0.at() + f_h.at());
-
-      set_bdry_conds(k1[0],k1[1]);
-      
-      cout << "k1 done!" << endl;
-
-      //debug Euler
-      u0 = k1[0]; v0 = k1[1]; h0 = k1[2]; continue;
-           
-      //k2
-      f_u.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
-      f_v.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
-      f_h.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
-
-      k2[0].interpolate(f_u.at());
-      k2[1].interpolate(f_v.at());
-      k2[2].interpolate(f_h.at());
- 
-      set_bdry_conds(k2[0],k2[1]);
+    rbf::conical::set_n(5);     
+    interpolator<RBF_TYPE> u(the_bvp);
+    u.precompute_ev();
+    interpolator<RBF_TYPE> u0 = u, u1 = u;
 
-      cout << "k2 done!" << endl;
-
-      //k3      
-      f_u.set_interps(u0+(k2[0]/2), v0+(k2[1]/2), h0+(k2[2]/2));
-      f_v.set_interps(u0+(k2[0]/2), v0+(k2[1]/2), h0+(k2[2]/2));
-      f_h.set_interps(u0+(k2[0]/2), v0+(k2[1]/2), h0+(k2[2]/2));
-
-      k3[0].interpolate(f_u.at());
-      k3[1].interpolate(f_v.at());
-      k3[2].interpolate(f_h.at());
-      
-      set_bdry_conds(k3[0],k3[1]);
-      
-      cout << "k3 done!" << endl;
-
-      //k4
-      f_u.set_interps(u0+k3[0], v0+k3[1], h0+k3[2]);
-      f_v.set_interps(u0+k3[0], v0+k3[1], h0+k3[2]);
-      f_h.set_interps(u0+k3[0], v0+k3[1], h0+k3[2]);
-
-      k4[0].interpolate(f_u.at());
-      k4[1].interpolate(f_v.at());
-      k4[2].interpolate(f_h.at());
-      
-      set_bdry_conds(k4[0],k4[1]);
-
-      cout << "k4 done!" << endl;
-     
-      //Grand finale
-      u0.interpolate((u0 + (k1[0] + 2*k2[0] + 2*k3[0] + k4[0])/6).at());
-      v0.interpolate((v0 + (k1[1] + 2*k2[1] + 2*k3[1] + k4[1])/6).at());
-      h0.interpolate((h0 + (k1[2] + 2*k2[2] + 2*k3[2] + k4[2])/6).at());
-
-      set_bdry_conds(u0,v0);
+    //Main loop
+    size_t maxiter = 1000;
+    for(size_t n = 1; n <= maxiter; n++)
+    {
+      u0 = u1;
+      u1 = u;
+      cout << "Iteration: " << n << endl;
+      u.set_f( (u0 - 2*u1).at() );
     }
 
-    return 0;
   }
   catch(error exc){
     utils::show_exception(exc);
@@ -249,133 +87,4 @@
 
 }
 
-//******************* Function definitions ***************************
 
-template<typename RBF>
-void Fgen<RBF>::set_interps(const interpolator<RBF>& u0, 
-			    const interpolator<RBF>& v0,
-			    const interpolator<RBF>& h0){
-  set_u(u0);
-  set_v(v0);
-  set_h(h0);
-}
-
-template<typename RBF>
-void Fgen<RBF>::set_u(const interpolator<RBF>& u0)
-{
-  u = u0;
-}
-
-template<typename RBF>
-void Fgen<RBF>::set_v(const interpolator<RBF>& v0)
-{
-  v = v0;
-}
-
-template<typename RBF>
-void Fgen<RBF>::set_h(const interpolator<RBF>& h0)
-{
-  h = h0;
-}
-
-template<typename RBF>
-void Fgen<RBF>::set_dt(double dt_in)
-{
-  dt = dt_in;
-}
-
-
-template<typename RBF>
-interp_values Fu<RBF>::at() const
-{
-  return   dt*(u()*u.d(1) +  v()*u.d(2) +  g*h.d(1));
-}
-
-template<typename RBF>
-interp_values Fv<RBF>::at() const
-{
-  return   dt*(u()*v.d(1) 
-	       +  v()*v.d(2) 
-	       +  g*h.d(2));
-}
-
-template<typename RBF>
-interp_values Fh<RBF>::at() const
-{
-  return   dt*(u()*h.d(1) + h()*u.d(1) +   v()*h.d(2) + h()*v.d(2));
-}
- 
-template<typename RBF>
-void save_interp(const shared_ptr<domain> Omega, 
-		 const interpolator<RBF>& u, 
-		 size_t n, string which_interp)
-{
-  using namespace std;
-  string filename = "results/" + which_interp;
-  stringstream ss;
-  ss << std::setw(5) << setfill('0') << n;
-  string n_string;
-  ss >> n_string;
-  filename += n_string + ".map";
-  ofstream ofs(filename.c_str());    
-  {
-    slice s(1,2);
-    matrix M(Omega -> get_interior().size() +
-	     Omega -> get_boundary().size(), 3);
-    size_t k = 1;
-    for(set<point>::const_iterator i = Omega -> get_interior().begin();
-	i != Omega -> get_interior().end(); i++){
-      M(k,s) = *i;
-      M(k,3) = u(*i);
-      k++;
-    }
-    for(set<point>::const_iterator i = Omega -> get_boundary().begin();
-	i != Omega -> get_boundary().end(); i++){
-      M(k,s) = *i;
-      M(k,3) = u(*i);
-      k++;
-    }
-    ofs << M;
-  }
-
-
-  if(false){ //debug
-    double r, th, N, M;
-    N = 30; M  = 70;
-    linalg::matrix out(static_cast<size_t>(N),static_cast<size_t>(M));
-    size_t i,j;
-    for(r = 0, j=1;  j <= N; r += 1/(N-1),j++){
-      for(th = 0,i=1; i <= M; th += 2*pi/(M-1),i++){
-	linalg::vector p(2); 
-	p(1) = r*cos(th); 
-	p(2) = r*sin(th);
-	out(j,i) = u(p);
-      }
-    }
-  
-    ofs << out;
-  }//debug
-}
-
-
-template<typename RBF>
-void set_bdry_conds(interpolator<RBF> & u, interpolator<RBF> & v)
-{
-  bvp::bdry_values u_bdry(u.at()), v_bdry(v.at());
-  bvp::normals N=u.get_normals(), M=v.get_normals();
-  if(N != M){
-    failure exc;
-    exc.reason="u and v have incompatible domains. This should never"
-      " happen.";
-    exc.line=__LINE__; exc.file=__FILE__; throw exc;
-  }
-
-  //Project the vector (u,v) onto the perp of (N(1),N(2)) = (-N(2),N(1)).
-  bvp::bdry_values u_bdry_new = u_bdry*N(2)*N(2) - v_bdry*N(1)*N(2);
-  bvp::bdry_values v_bdry_new = v_bdry*N(1)*N(1) - u_bdry*N(1)*N(2);
-
-  u.set_bdry_values(u_bdry_new);
-  v.set_bdry_values(v_bdry_new);
-}
-
-

new file mode 100644
--- /dev/null
+++ b/results/showvecs.m
@@ -0,0 +1,14 @@
+N = 21810;
+skip = 10;
+
+for i = skip:skip:N
+  nameu = "u"; namev = "v";
+  istr = sprintf("%05d",i);
+  filenameu = strcat(nameu,istr,".map");
+  filenamev = strcat(namev,istr,".map");
+  u = load(filenameu);
+  v = load(filenamev);
+  x = u(:,1); y = u(:,2);
+  u = u(:,3); v = v(:,3);
+  quiver(x,y,u,v,0.15);
+endfor
\ No newline at end of file