--- a/main.cpp
+++ b/main.cpp
@@ -36,7 +36,7 @@
 const double g = 9.8; // m/s^2
 
 template<typename RBF>
-class Fgen : public realfunc{
+class Fgen {
 public:
   void set_interps(const interpolator<RBF>& u0, 
 		   const interpolator<RBF>& v0,
@@ -59,7 +59,7 @@
   using Fgen<RBF>::h;
   using Fgen<RBF>::dt;
 public:
-  double at(const point& p) const;
+  bvp::interp_values at() const;
 };
 
 template<typename RBF>
@@ -69,7 +69,7 @@
   using Fgen<RBF>::h;
   using Fgen<RBF>::dt;
 public:
-  double at(const point& p) const;
+  bvp::interp_values at() const;
 };
 
 template<typename RBF>
@@ -79,9 +79,10 @@
   using Fgen<RBF>::h;
   using Fgen<RBF>::dt;
 public:
-  double at(const point& p) const;
+  bvp::interp_values at() const;
 };
 
+//FIXME: This iteration must work on interpolators without boundary
 class iter_neumann : public bdry_diff_op{
 public:
   //If this is u's, set 1, for v, set 2.
@@ -164,15 +165,17 @@
       h1, h0(Omega, h_init);
 
     //Precompute some RBFs
-    u0.precompute_at(); v0.precompute_at(), h0.precompute_at();
-    u0.precompute_d(); v0.precompute_d(), h0.precompute_d();
+    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> >
-      k1(3,h0), k2(3,h0), k3(3,h0), k4(3,h0);
+      // 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); 
 
     Fu<RBF_TYPE> f_u;
     Fv<RBF_TYPE> f_v;
@@ -182,8 +185,6 @@
     f_v.set_dt(dt);
     f_h.set_dt(dt);
 
-
-
     //main loop
     size_t maxiter = 1000;
     for(size_t i = 1; i <= maxiter; i++){ 
@@ -199,9 +200,6 @@
       f_v.set_interps(u0,v0,h0);
       f_h.set_interps(u0,v0,h0);
       
-      k1[0].set_f(f_u);
-      k1[1].set_f(f_v);
-      k1[2].set_f(f_h);
       cout << "k1 done!" << endl;
 
       bdry_iter(k1[0],k1[1],Omega);
@@ -211,10 +209,6 @@
       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].set_f(f_u);
-      k2[1].set_f(f_v);
-      k2[2].set_f(f_h);
-
       bdry_iter(k2[0],k2[1],Omega);
 
       cout << "k2 done!" << endl;
@@ -224,10 +218,6 @@
       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].set_f(f_u);
-      k3[1].set_f(f_v);
-      k3[2].set_f(f_h);
-
       bdry_iter(k3[0],k3[1],Omega);
       
       cout << "k3 done!" << endl;
@@ -237,10 +227,6 @@
       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].set_f(f_u);
-      k4[1].set_f(f_v);
-      k4[2].set_f(f_h);
-
       bdry_iter(k4[0],k4[1],Omega);
 
       cout << "k4 done!" << endl;
@@ -262,7 +248,7 @@
   }
   catch(error exc){
     utils::show_exception(exc);
-    return 1;
+    return ;
   }
 
 }
@@ -304,21 +290,21 @@
 
 
 template<typename RBF>
-double Fu<RBF>::at(const point& p) const
+interp_values Fu<RBF>::at() const
 {
-  return   dt*(u(p)*u.d(p,1) +  v(p)*u.d(p,2) +  g*h.d(p,1));
+  return   dt*(u()*u.d(1) +  v()*u.d(2) +  g*h.d(1));
 }
 
 template<typename RBF>
-double Fv<RBF>::at(const point& p) const
+interp_values Fv<RBF>::at() const
 {
-  return   dt*(u(p)*v.d(p,1) +  v(p)*v.d(p,2) +  g*h.d(p,2));
+  return   dt*(u()*v.d(1) +  v()*v.d(2) +  g*h.d(2));
 }
 
 template<typename RBF>
-double Fh<RBF>::at(const point& p) const
+interp_values Fh<RBF>::at() const
 {
-  return   dt*(u(p)*h.d(p,1) + h(p)*u.d(p,1) +   v(p)*h.d(p,2) + h(p)*v.d(p,2));
+  return   dt*(u()*h.d(1) + h()*u.d(1) +   v()*h.d(2) + h()*v.d(2));
 }
 
 iter_neumann::iter_neumann(size_t u_or_v_in)