diff --git a/poisson/double_buffer.cc b/poisson/double_buffer.cc
index e6ecf0f..ac5be0d 100644
--- a/poisson/double_buffer.cc
+++ b/poisson/double_buffer.cc
@@ -1,19 +1,25 @@
 /* -------------------------------------------------------------------------- */
 #include "double_buffer.hh"
 #include "grid.hh"
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 DoubleBuffer::DoubleBuffer(int m, int n)
     : m_current(new Grid(m, n)), m_old(new Grid(m, n)) {}
 
 /* -------------------------------------------------------------------------- */
 Grid & DoubleBuffer::current() { return *m_current; }
 
 /* -------------------------------------------------------------------------- */
 Grid & DoubleBuffer::old() { return *m_old; }
 
+/* -------------------------------------------------------------------------- */
+void DoubleBuffer::resize(int m, int n) {
+  m_current->resize(m, n);
+  m_old->resize(m, n);
+}
+
 /* -------------------------------------------------------------------------- */
 void DoubleBuffer::swap() {
   m_current.swap(m_old);
 }
diff --git a/poisson/double_buffer.hh b/poisson/double_buffer.hh
index 67a1405..88efba9 100644
--- a/poisson/double_buffer.hh
+++ b/poisson/double_buffer.hh
@@ -1,24 +1,27 @@
 #ifndef DOUBLE_BUFFER
 #define DOUBLE_BUFFER
 
 /* -------------------------------------------------------------------------- */
 #include <memory>
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 /* -------------------------------------------------------------------------- */
 
 class DoubleBuffer {
 public:
-  DoubleBuffer(int m, int n);
+  DoubleBuffer(int m = 0, int n = 0);
 
   Grid & current();
   Grid & old();
 
+  // resize the grid
+  void resize(int m, int n);
+
   void swap();
 private:
 
   std::unique_ptr<Grid> m_current;
   std::unique_ptr<Grid> m_old;
 };
 
 #endif /* DOUBLE_BUFFER */
diff --git a/poisson/grid.cc b/poisson/grid.cc
index 4e9fb84..158e7ce 100644
--- a/poisson/grid.cc
+++ b/poisson/grid.cc
@@ -1,15 +1,22 @@
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 /* -------------------------------------------------------------------------- */
 #include <algorithm>
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 Grid::Grid(int m, int n) : m_m(m), m_n(n), m_storage(m * n) { clear(); }
 
 /* -------------------------------------------------------------------------- */
 void Grid::clear() { std::fill(m_storage.begin(), m_storage.end(), 0.); }
 
+/* -------------------------------------------------------------------------- */
+void Grid::resize(int m, int n) {
+  m_m = m;
+  m_n = n;
+  m_storage.resize(m*n);
+}
+
 /* -------------------------------------------------------------------------- */
 int Grid::m() const { return m_m; }
 int Grid::n() const { return m_n; }
diff --git a/poisson/grid.hh b/poisson/grid.hh
index ad872cb..23c53e1 100644
--- a/poisson/grid.hh
+++ b/poisson/grid.hh
@@ -1,29 +1,32 @@
 #ifndef GRID_HH
 #define GRID_HH
 
 /* -------------------------------------------------------------------------- */
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 class Grid {
 public:
-  Grid(int m, int n);
+  Grid(int m = 0, int n = 0);
 
   /// access the value [i][j] of the grid
   inline float & operator()(int i, int j) { return m_storage[i * m_n + j]; }
   inline const float & operator()(int i, int j) const {
     return m_storage[i * m_n + j];
   }
 
+  // resize the grid
+  void resize(int m, int n);
+
   /// set the grid to 0
   void clear();
 
   int m() const;
   int n() const;
 
 private:
   int m_m, m_n;
   std::vector<float> m_storage;
 };
 
 #endif /* GRID_HH */
diff --git a/poisson/simulation.cc b/poisson/simulation.cc
index f954d4e..293cda4 100644
--- a/poisson/simulation.cc
+++ b/poisson/simulation.cc
@@ -1,65 +1,65 @@
 /* -------------------------------------------------------------------------- */
 #include "simulation.hh"
 /* -------------------------------------------------------------------------- */
 #include <cmath>
 #include <iostream>
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 Simulation::Simulation(int m, int n)
     : m_global_m(m), m_global_n(n), m_epsilon(1e-7), m_h_m(1. / m),
       m_h_n(1. / n), m_grids(m, n), m_f(m, n),
       m_dumper(new DumperASCII(m_grids.old())) {}
 
 /* -------------------------------------------------------------------------- */
 void Simulation::set_initial_conditions() {
   for (int i = 0; i < m_global_m; i++) {
     for (int j = 0; j < m_global_n; j++) {
       m_f(i, j) = -2. * 100. * M_PI * M_PI * std::sin(10. * M_PI * i * m_h_m) *
                   std::sin(10. * M_PI * j * m_h_n);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 std::tuple<float, int> Simulation::compute() {
   int s = 0;
   float l2 = 0;
   do {
     l2 = compute_step();
 
     m_grids.swap();
 
     m_dumper->dump(s);
     ++s;
   } while (l2 > m_epsilon);
 
   return std::make_tuple(std::sqrt(l2), s);
 }
 
 /* -------------------------------------------------------------------------- */
 void Simulation::set_epsilon(float epsilon) { m_epsilon = epsilon * epsilon; }
 
 /* -------------------------------------------------------------------------- */
 float Simulation::epsilon() const { return std::sqrt(m_epsilon); }
 
 /* -------------------------------------------------------------------------- */
 float Simulation::compute_step() {
   float l2 = 0.;
 
   Grid & u = m_grids.current();
   Grid & uo = m_grids.old();
 
   for (int i = 1; i < m_global_m - 1; i++) {
-    for (int j = 1; j < m_global_m - 1; j++) {
+    for (int j = 1; j < m_global_n - 1; j++) {
       // computation of the new step
       u(i, j) = 0.25 * (uo(i - 1, j) + uo(i + 1, j) + uo(i, j - 1) +
                         uo(i, j + 1) - m_f(i, j) * m_h_m * m_h_n);
 
       // L2 norm
       l2 += (uo(i, j) - u(i, j)) * (uo(i, j) - u(i, j));
     }
   }
 
   return l2;
 }