diff --git a/src/action.hpp b/src/action.hpp
index feea85e3e4c98f8686b99aed20033c675047eae2..70831e3ffc2993bb54567c73036847a079e811ab 100644
--- a/src/action.hpp
+++ b/src/action.hpp
@@ -8,9 +8,36 @@
 template <typename Scalar, typename Shape>
 inline Scalar S_2(Array<Scalar, Shape> Q)
 {
-    Q(0u,0u) = 1;
-    Q(0u,1u) = 2;
-    Q(1u,1u) = 3;
+    size_t nx = Q.dim(0);
+    size_t ny = Q.dim(1);
+
+    for (size_t i = 0; i < nx; ++i) {
+        for (size_t j = 0; j < ny; ++j) {
+            auto Q1 = Q.part(i,j).matrix();
+
+            /* Sum over neighbors */
+            for (int di = -1; di <= 1; ++di) {
+                for (int dj = -1; dj <= 1; ++dj) {
+                    if (i == 0 && di < 0)
+                        continue;
+                    if (i == nx-1 && di > 0)
+                        continue;
+                    if (j == 0 && dj < 0)
+                        continue;
+                    if (j == ny-1 && dj > 0)
+                        continue;
+                    if (di == 0 && dj == 0)
+                        continue;
+
+                    size_t i2 = i + di;
+                    size_t j2 = j + dj;
+
+                    auto Q2 = Q.part(i2,j2).matrix();
+                }
+            }
+        }
+    }
+
     return Scalar(0);
 }
 
diff --git a/src/array.hpp b/src/array.hpp
index 1ed483fc2685ff614a334862ce230b8b69ac0627..b06a7fbd716ce5c99573a470eb5f4780b0676c98 100644
--- a/src/array.hpp
+++ b/src/array.hpp
@@ -21,6 +21,7 @@ struct Shape<>
     static const size_t ndim = 0;
     static const size_t size = 1;
     static const bool fixed = true;
+    typedef Shape<> Tail;
     static constexpr size_t dim(size_t i) { return 0; }
     static constexpr size_t stride(size_t i) { return 1; }
 };
@@ -192,7 +193,8 @@ private:
     template <size_t axis>
     static constexpr Eigen::Index eigen_shape()
         {
-            size_t n = (axis == 0) ? Shape::head : Shape::Tail::head;
+            /* bad Shape handled in matrix() */
+            size_t n = (axis == 0 || Shape::ndim < 2) ? Shape::head : Shape::Tail::head;
             return (n == Dynamic) ? Eigen::Dynamic : n;
         }
 
@@ -289,7 +291,7 @@ public:
         {
             constexpr size_t nidxs = sizeof...(Idx);
             static_assert(nidxs < Shape::ndim,
-                          "number of indices must be small than the number of dimensions");
+                          "number of indices must be smaller than the number of dimensions");
 
             size_t offset = index(idxs...);
             return {data_ + offset, Base::size() - offset};
diff --git a/src/main.cpp b/src/main.cpp
index 199216dea96bca6e746fbe34c8177add0849b2bc..ee71d97e049e5fd73aea4d963523673307e5430d 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -24,11 +24,11 @@ void dump(Array<Scalar,Shape> Q)
 
 int main()
 {
-    std::vector<double> v(1*2*2*2);
-    Array<double,Shape<Dynamic,Dynamic,2,2> > Q0(v.data(), v.size(), {1,2,2,2});
+    std::vector<double> v(10*10*2*2);
+    Array<double,Shape<Dynamic,Dynamic,2,2> > Q0(v.data(), v.size(), {10,10,2,2});
     auto Q = Q0.part(0u,1u);
 
-    S_2(Q);
+    S_2(Q0);
 
     std::cout << "size = " << Q.size() << std::endl;