--- a/main.cpp
+++ b/main.cpp
@@ -28,7 +28,7 @@
   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);
+    return dt2*L(f,p) - f(p);
   }
 private:
   Laplacian L;
@@ -40,6 +40,9 @@
   double at(const point&) const{return 0;};
 };
 
+void save_iteration(const interpolator<RBF_TYPE>& u,
+                    size_t n);
+
 int main()
 {  
   gsl_set_error_handler(&error_handling::errorHandler);
@@ -52,8 +55,7 @@
       
     shared_ptr<domain> Omega(new domain("data/circ_intr.matrix",
                                         "data/circ_bdry.matrix",
-                                        "data/circ_nrml.matrix"));
-    
+                                        "data/circ_nrml.matrix"));    
     zero_func g;
 
     //timestep
@@ -70,9 +72,10 @@
     interpolator<RBF_TYPE> u0 = u, u1 = u;
 
     //Main loop
-    size_t maxiter = 1000;
+    size_t maxiter = 500;
     for(size_t n = 1; n <= maxiter; n++)
     {
+      save_iteration(u,n);
       u0 = u1;
       u1 = u;
       cout << "Iteration: " << n << endl;
@@ -87,4 +90,16 @@
 
 }
 
-
+void 
+save_iteration(const interpolator<RBF_TYPE>& u,
+               size_t n)
+{
+  using namespace std;
+  stringstream ss;
+  string filename = "results/u";
+  ss << filename << setw(5) << setfill('0') << n << ".map";
+  ss >> filename;
+  ofstream ofs(filename.c_str());
+  ofs << u;
+  ofs.close();
+}