--- a/Makefile
+++ b/Makefile
@@ -39,4 +39,8 @@
 main.o: main.cpp $(HEADERFILES)
 
 clean:
-	rm -f *.o foo
\ No newline at end of file
+	rm -f *.o foo
+
+cleandata:
+	rm -f data/*.matrix data/*.map data/*.vector \
+	 	results/*.matrix results/*.map results/*.png results/*.vector
\ No newline at end of file

--- a/include/interpolator.hpp
+++ b/include/interpolator.hpp
@@ -65,11 +65,27 @@
     double d2(const point &p, size_t k1, size_t k2) const;
     //@}
 
+    /** @name Precomputation of RBFs
+     *
+     * The following functions will precompute upon request all the
+     * interpolator's RBFs at all points of the domain. This may speed
+     * up some computations at the expense of increasing memory
+     * usage.
+     */
+    //@{
+    ///Precompute all the function evaluation of the RBFs.
+    void precompute_at();
+    ///Precompute all relevant first derivatives.
+    void precompute_d() ;
+    ///Precompute all relevant second derivatives.
+    void precompute_d2();
+    //@}
+
     /** @name Partial redefinitions
      *
      * These functions allow for partial redefinition of the BVP as
      * equired for the additive Schwartz domain decomposition method,
-     * a nd for other methods. They do not factor a matrix again.
+     * and for other methods. They do not factor a matrix again.
      */
     //@{
     void set_f(const realfunc &f);
@@ -81,7 +97,8 @@
     /** @name Linear arithmetic operators
      *
      *  These functions return a new interpolator. They are pointwise
-     *  linear operations.
+     *  linear operations and work by manipulating interpolator
+     *  coefficients.
      */
     //@{
     /// Needs two operators on the same domain.
@@ -94,36 +111,35 @@
     interpolator<RBF> operator/(double a) const;
     //@}
     
-    //debug
-    mutable map<std::vector<size_t>, matrix> precomp_rbfs;
   private:
-    //Once the matrix is defined, this function inverts it.
+    ///Once the matrix is defined, this function inverts it.
     void computecoeffs(); 
 	 
-    //Perform the actual interpolation.
+    ///Perform the actual interpolation.
     void init(shared_ptr<linear_BVP2> bvp);
 
     shared_ptr<linear_BVP2> thebvp;
 
-    //Number of interior points. 
+    ///Number of interior points. 
     size_t n; 
-    //Number of boundary points.
+    ///Number of boundary points.
     size_t m; 
 
-    //The matrix to invert.
+    ///The matrix to invert.
     matrix M; 
 
-    //Is the interpolator ready for use?
+    ///Is the interpolator ready for use?
     bool initted;
     void not_initted(int line, string file) const; //Exception thrower
     
-    //Coefficients of the RBFs
+    ///Coefficients of the RBFs
     linalg::vector coeffs;
-    //The RBFS
+    ///The RBFS
     std::vector<RBF> rbfs;
-    //Their hash
+    ///Their hash
     size_t rbfs_hash;    
 
+    ///An RBF hasher
     size_t hash_value(const std::vector<RBF>& rbfs_in);
        
     /** \short Precomputed RBFs 
@@ -139,10 +155,17 @@
      * 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?
+    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;
+    mutable map<std::vector<size_t>, map<point, double> >
+    precomp_values;
+    /// Invert phis\coeffs and return result vector in a map.
     map<point, double> precompute_values(const matrix& phis) const;
   };
 

--- a/include/linalg.hpp
+++ b/include/linalg.hpp
@@ -266,7 +266,7 @@
     vector operator+(const vector& w) const;
     ///Vector subtration.
     vector operator-(const vector& w) const;
-    ///xDot product.
+    ///Dot product.
     double operator*(const vector& w) const; 
     ///Computes vM where v is treated as a row vector.
     vector operator*(const matrix& M) const; 

--- a/interpolator.cpp
+++ b/interpolator.cpp
@@ -12,9 +12,6 @@
 #include <boost/functional/hash.hpp>
 #include <boost/shared_ptr.hpp>
 
-//debug
-#include <iostream>
-
 namespace bvp{ 
   using boost::shared_ptr;
 
@@ -191,8 +188,14 @@
     }
    
     computecoeffs();
+    
+    ev_precomp = false;
+    d1_precomp = false;
+    d2_precomp = false;
+
+    rbfs_hash = hash_value(rbfs);
+
     initted = true;
-    rbfs_hash = hash_value(rbfs);
   }
 
   template<typename RBF>
@@ -234,17 +237,17 @@
     }
     return out;
   }
-
+  
   template<typename RBF>
-  double interpolator<RBF>::at(const point& p) const{
+  void interpolator<RBF>::precompute_at(){
     if(!initted){
       not_initted(__LINE__, __FILE__);
     }
-    std::vector<size_t> alpha; //empty vector
-    
-    //Have the RBFs been precomputed yet?
-    if(precomp_rbfs.find(alpha) == precomp_rbfs.end()){
-      //No, so precompute the RBFs at all points of the domain
+
+    if(!ev_precomp){
+      std::vector<size_t> alpha; //empty vector
+      
+      //Precompute the RBFs at all points of the domain
       matrix phis(n+m,n+m);
       set<point>::iterator I;
       size_t i;
@@ -262,19 +265,65 @@
       precomp_values[alpha] = precompute_values(phis);
     }
 
-    
+    ev_precomp = true;
+  }
+
+  template<typename RBF>
+  double interpolator<RBF>::at(const point& p) const{
+    if(!initted){
+      not_initted(__LINE__, __FILE__);
+    }
+    std::vector<size_t> alpha; //empty vector
+
     //Are we evaluating at a precomputed point in the domain?
-    if(precomp_values[alpha].find(p) != precomp_values[alpha].end())
+    if(ev_precomp and 
+       precomp_values[alpha].find(p) != precomp_values[alpha].end())
       return precomp_values[alpha][p];
+     
     
     //Else, must compute the value
     double result = 0;
     for(size_t i = 1; i <= coeffs.size(); i++)
       result += coeffs(i)*rbfs[i-1].at(p);
-    
-    
+        
     return result;
   }
+  
+  template<typename RBF>
+  void interpolator<RBF>::precompute_d()
+  {
+    if(!initted){
+      not_initted(__LINE__, __FILE__);
+    }
+
+    //Have the RBFs been precomputed yet?
+    if(!d1_precomp){
+      //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);
+      
+	  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;
+	precomp_values[alpha] = precompute_values(phis);
+      }
+    }
+    d1_precomp = true;
+  }
 
   template<typename RBF>
   double interpolator<RBF>::d(const point& p, size_t k) const{
@@ -282,31 +331,10 @@
       not_initted(__LINE__, __FILE__);
     }
     std::vector<size_t> alpha(k); alpha[k-1]++;
-
-    //Have the RBFs been precomputed yet?
-    if( precomp_rbfs.find(alpha) == precomp_rbfs.end()){
-      //No, so precompute the RBFs at all points of the domain
-      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);
-      }
-
-      precomp_rbfs[alpha] = phis;
-      precomp_values[alpha] = precompute_values(phis);
-    }
     
     //Are we evaluating at a precomputed point in the domain?
-    if(precomp_values[alpha].find(p) != precomp_values[alpha].end())
+    if(d1_precomp and 
+       precomp_values[alpha].find(p) != precomp_values[alpha].end())
       return precomp_values[alpha][p];
     
     //Else, must compute the value
@@ -319,34 +347,48 @@
   }
 
   template<typename RBF>
+  void interpolator<RBF>::precompute_d2()
+  {
+    if(!initted)
+      not_initted(__LINE__, __FILE__);
+    if(!d2_precomp){
+      size_t dim = thebvp -> get_domain() -> get_dimension();
+      //precompute all second derivatives of t
+      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(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;
+	    precomp_values[alpha] = precompute_values(phis);
+	  
+	}
+      }
+    }
+    d2_precomp = true;
+  }
+
+  template<typename RBF>
   double interpolator<RBF>::d2(const point& p, size_t k1, size_t k2) const{
     if(!initted){
       not_initted(__LINE__, __FILE__);
     }
     std::vector<size_t> alpha(k1>k2?k1:k2); alpha[k1-1]++; alpha[k2-1]++;
-
-    //Have the RBFs been precomputed yet?
-    if( precomp_rbfs.find(alpha) == precomp_rbfs.end()){
-      //No, so precompute the RBFs at all points of the domain
-      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);
-      
-      precomp_rbfs[alpha] = phis;
-      precomp_values[alpha] = precompute_values(phis);
-    }
     
     //Are we evaluating at a precomputed point in the domain?
-    if(precomp_values[alpha].find(p) != precomp_values[alpha].end())
+    if(d2_precomp and 
+       precomp_values[alpha].find(p) != precomp_values[alpha].end())
       return precomp_values[alpha][p];
     
     //Else, must compute the value

--- a/main.cpp
+++ b/main.cpp
@@ -163,11 +163,10 @@
       v1, v0(v_bvp_init), 
       h1, h0(Omega, h_init);
 
-    //FIXME: put proper precomp function here.
-    point p = *(Omega -> get_interior().begin());
-    u0(p); v0(p); h0(p); 
-    u0.d(p,1); v0.d(p,1); h0.d(p,1); 
-    u0.d(p,2); v0.d(p,2); h0.d(p,2); 
+    //Precompute some RBFs
+    u0.precompute_at(); v0.precompute_at(), h0.precompute_at();
+    u0.precompute_d(); v0.precompute_d(), h0.precompute_d();
+
 
     u1 = u0; v1 = v0; h1 = h0;
 
@@ -186,7 +185,7 @@
 
 
     //main loop
-    size_t maxiter = 10;
+    size_t maxiter = 1000;
     for(size_t i = 1; i <= maxiter; i++){ 
       cout << "Now on iteration #" << i << endl;
       if(i % 1 == 0){
@@ -203,6 +202,7 @@
       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);
            
@@ -217,6 +217,8 @@
 
       bdry_iter(k2[0],k2[1],Omega);
 
+      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));
@@ -227,6 +229,8 @@
       k3[2].set_f(f_h);
 
       bdry_iter(k3[0],k3[1],Omega);
+      
+      cout << "k3 done!" << endl;
 
       //k4
       f_u.set_interps(u0+k3[0], v0+k3[1], h0+k3[2]);
@@ -238,6 +242,8 @@
       k4[2].set_f(f_h);
 
       bdry_iter(k4[0],k4[1],Omega);
+
+      cout << "k4 done!" << endl;
      
       //Grand finale
       u1.set_f(u0 + (k1[0] + 2*k2[0] + 2*k3[0] + k4[0])/6);