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

--- a/data/sw_eqs_init.m
+++ b/data/sw_eqs_init.m
@@ -1,5 +1,5 @@
-n = 12
-m = 10;
+n = 15
+m = 12;
 r = linspace(0,1,n);
 
 rrintr = 0;

--- a/include/interpolator.hpp
+++ b/include/interpolator.hpp
@@ -190,6 +190,10 @@
     ///Perform the actual interpolation.
     void init(shared_ptr<linear_BVP2> bvp);
 
+    ///If the coefficients change, need to recompute precomputed values 
+    void recompute_values_vec();
+    
+    ///BVP associated to this interpolator
     shared_ptr<linear_BVP2> thebvp;
 
     ///Number of interior points. 

--- a/interpolator.cpp
+++ b/interpolator.cpp
@@ -156,15 +156,7 @@
     coeffs = M.inv(values.v);
    
     //Precompute the relevant values
-    {
-      typename 
-	map<std::vector<size_t>, matrix>::const_iterator I;
-      for(I = precomp_rbfs.begin(); I != precomp_rbfs.end(); I++)
-	precomp_values_vec[I -> first] = (I -> second)*coeffs;
-      
-      //Pointwise map evaluations we will not compute again.
-      precomp_values.clear();
-    }
+    recompute_values_vec();
   }
 
   template<typename RBF>
@@ -486,6 +478,7 @@
     
     interpolator<RBF> out = *this;
     out.coeffs = (this -> coeffs) + u.coeffs;
+    out.recompute_values_vec();
        
     return out;
   }
@@ -505,6 +498,7 @@
 
     interpolator<RBF> out = *this;
     out.coeffs = (this -> coeffs) - u.coeffs;
+    out.recompute_values_vec();
        
     return out;
   }
@@ -514,6 +508,7 @@
   {
     interpolator<RBF> u = *this;
     u.coeffs = (this -> coeffs)*a;
+    u.recompute_values_vec();
     return u;
   }
 
@@ -522,6 +517,7 @@
   {
     interpolator<RBF> u = *this;
     u.coeffs = (this -> coeffs)*(1/a);
+    u.recompute_values_vec();
     return u;
   }
 
@@ -587,16 +583,19 @@
     slice s1(1,n), s2(n+1,n+m);
     rhs(s2) = b_new.v;
     coeffs = precomp_rbfs[alpha].inv(rhs);
-    //FIXME: Don't repeat this code
-    {
-      typename 
-	map<std::vector<size_t>, matrix>::const_iterator I;
-      for(I = precomp_rbfs.begin(); I != precomp_rbfs.end(); I++)
-	precomp_values_vec[I -> first] = (I -> second)*coeffs;
+    
+    recompute_values_vec();
+  }
+
+  template<typename RBF>
+  void interpolator<RBF>::recompute_values_vec(){
+    typename 
+      map<std::vector<size_t>, matrix>::const_iterator I;
+    for(I = precomp_rbfs.begin(); I != precomp_rbfs.end(); I++)
+      precomp_values_vec[I -> first] = (I -> second)*coeffs;
       
-      //Pointwise map evaluations we will not compute again.
-      precomp_values.clear();
-    }
+    //Pointwise map evaluations we will not compute again.
+    precomp_values.clear();
   }
 
   template<typename RBF>

--- a/main.cpp
+++ b/main.cpp
@@ -18,7 +18,7 @@
 using namespace linalg;
 using namespace bvp;
 
-#define RBF_TYPE rbf::conical
+#define RBF_TYPE rbf::thin_plate_spline
 
 //Will solve the system
 //
@@ -122,28 +122,23 @@
       u_bvp_init(new linear_BVP2(Omega, I, D, Z, Z)),
       v_bvp_init(new linear_BVP2(Omega, I, D, Z, Z));
       
-   
-    double dt = 1;
-
+    
     //init the interps.
     using namespace rbf;
     conical::set_n(5); thin_plate_spline::set_n(6); 
-    c_infty_rbf::set_epsilon(0.01);
+    c_infty_rbf::set_epsilon(10);
 
     cout << "Initialising... please wait." << endl;
 
     interpolator<RBF_TYPE> 
-      u1, u0(u_bvp_init), 
-      v1, v0(v_bvp_init), 
-      h1, h0(Omega, h_init);
+      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();
 
-
-    u1 = u0; v1 = v0; h1 = h0;
-
     //Intermediate interpolators for RK4
     std::vector<interpolator<RBF_TYPE> >
       // FIXME: make the copy ctor for interps make different copies
@@ -154,6 +149,9 @@
     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);
@@ -177,7 +175,7 @@
       k1[1].interpolate(f_v.at());
       k1[2].interpolate(f_h.at());
 
-      //set_bdry_conds(k1[0],k1[1]);
+      set_bdry_conds(k1[0],k1[1]);
       
       cout << "k1 done!" << endl;
            
@@ -190,7 +188,7 @@
       k2[1].interpolate(f_v.at());
       k2[2].interpolate(f_h.at());
  
-      //set_bdry_conds(k2[0],k2[1]);
+      set_bdry_conds(k2[0],k2[1]);
 
       cout << "k2 done!" << endl;
 
@@ -203,7 +201,7 @@
       k3[1].interpolate(f_v.at());
       k3[2].interpolate(f_h.at());
       
-      //set_bdry_conds(k3[0],k3[1]);
+      set_bdry_conds(k3[0],k3[1]);
       
       cout << "k3 done!" << endl;
 
@@ -216,20 +214,16 @@
       k4[1].interpolate(f_v.at());
       k4[2].interpolate(f_h.at());
       
-      //set_bdry_conds(k4[0],k4[1]);
+      set_bdry_conds(k4[0],k4[1]);
 
       cout << "k4 done!" << endl;
      
       //Grand finale
-      u1.interpolate((u0 + (k1[0] + 2*k2[0] + 2*k3[0] + k4[0])/6).at());
-      v1.interpolate((v0 + (k1[1] + 2*k2[1] + 2*k3[1] + k4[1])/6).at());
-      h1.interpolate((h0 + (k1[2] + 2*k2[2] + 2*k3[2] + k4[2])/6).at());
-    
-      //set_bdry_conds(u1,v1);
-      
-      u0 = u1;
-      v0 = v1;
-      h0 = h1;
+      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);
     }
 
     return 0;
@@ -286,7 +280,9 @@
 template<typename RBF>
 interp_values Fv<RBF>::at() const
 {
-  return   dt*(u()*v.d(1) +  v()*v.d(2) +  g*h.d(2));
+  return   dt*(u()*v.d(1) 
+	       +  v()*v.d(2) 
+	       +  g*h.d(2));
 }
 
 template<typename RBF>