--- a/main.cpp
+++ b/main.cpp
@@ -1,4 +1,5 @@
 #include <iostream>
+#include <iomanip>
 #include <fstream>
 #include <map>
 #include <set>
@@ -18,7 +19,7 @@
 using namespace linalg;
 using namespace bvp;
 
-#define RBF_TYPE rbf::thin_plate_spline
+#define RBF_TYPE rbf::conical
 
 //Will solve the system
 //
@@ -35,6 +36,8 @@
 
 const double g = 9.8; // m/s^2
 
+const double pi = 3.141592653589793238462643383279502;
+
 template<typename RBF>
 class Fgen {
 public:
@@ -125,8 +128,8 @@
     
     //init the interps.
     using namespace rbf;
-    conical::set_n(5); thin_plate_spline::set_n(6); 
-    c_infty_rbf::set_epsilon(10);
+    conical::set_n(7); thin_plate_spline::set_n(6); 
+    c_infty_rbf::set_epsilon(0.01);
 
     cout << "Initialising... please wait." << endl;
 
@@ -143,27 +146,30 @@
     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;
     Fh<RBF_TYPE> f_h;
 
     //Negative because of the way the F's below are written.
-    double dt = -0.01;
+    double dt = -0.0001;
 
     f_u.set_dt(dt);
     f_v.set_dt(dt);
     f_h.set_dt(dt);
 
     //main loop, timestepping with RK4.
-    size_t maxiter = 1000;
+    size_t maxiter = 5000;
     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");
+ 	save_interp(Omega,u0,i,"u");
+ 	save_interp(Omega,v0,i,"v");
+ 	save_interp(Omega,h0,i,"h");
       }
 
       //k1 
@@ -171,13 +177,21 @@
       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());
+//       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));
@@ -298,38 +312,49 @@
 {
   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() +
+	     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;
+  }
 
-  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
 }