diff --git a/.ycm_extra_conf.py b/.ycm_extra_conf.py
index 156103100b8885036931a2cbb148adc8371a05a3..09a57a749c8438dc4200eeea34b9cb9bc1845af4 100644
--- a/.ycm_extra_conf.py
+++ b/.ycm_extra_conf.py
@@ -40,13 +40,9 @@ flags = [
     '-x',
     'c++',
     '-I',
-    '.',
-    '-isystem',
-    '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1',
-    '-isystem',
-    '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1/bits',
-    '-I',
-    'src',
+    './src',
+#    '-isystem',
+#    '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1',
 #    '-I',
 #    '/usr/include/c++/4.9.2',
 #    '-isystem',
diff --git a/src/fvm.hpp b/src/fvm.hpp
index 0b4cb1d09a63df689d7aa2c88908e7cff079c05d..db19aaae308840f45445a09be20500bf0138b346 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -24,6 +24,8 @@ class fvm_cell {
     /// the integral index type
     using size_type  = I;
 
+    using matrix_type = matrix<value_type, size_type>;
+
     /// the container used for indexes
     using index_type = memory::HostVector<size_type>;
     /// view into index container
@@ -40,10 +42,15 @@ class fvm_cell {
     /// build the matrix for a given time step
     void setup_matrx(value_type dt);
 
+    matrix_type& matrix()
+    {
+        return matrix_;
+    }
+
     private:
 
     /// the linear system for implicit time stepping of cell state
-    matrix<value_type, size_type> matrix_;
+    matrix_type matrix_;
 
     /// cv_areas_[i] is the surface area of CV i
     vector_type cv_areas_;
@@ -88,7 +95,6 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
                 );
             }
             cv_areas_[0] += math::area_sphere(soma->radius());
-            // d[0] += cv_areas_[0];
         }
         else if(auto cable = s->as_cable()) {
             // loop over each compartment in the cable
@@ -119,17 +125,10 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
                 auto ar = math::area_frustrum(halflen, right(c.radius), radius_center);
                 cv_areas_[j] += al;
                 cv_areas_[i] += ar;
-
-                // d[j] += al;
-                // d[i] += ar;
-                // l[i] -= alpha_ij;
-                // u[i] -= alpha_ij;
             }
         }
         else {
-            throw std::domain_error(
-                    "FVM lowering encountered unsuported segment type"
-            );
+            throw std::domain_error("FVM lowering encountered unsuported segment type");
         }
         ++seg_idx;
     }
@@ -158,19 +157,20 @@ void fvm_cell<T, I>::setup_matrx(T dt)
     //        .     .  .
     //       l[i] . . d[i]
     //
-    d(all) = 0;
-    for(auto i : d.range()) {
+
+    //d(all) = cv_areas_ + dt*(alpha_ + alpha_(p));
+    //d[0]  = cv_areas_[0];
+
+    d(all) = cv_areas_;
+    for(auto i=1; i<d.size(); ++i) {
         auto a = dt * alpha_[i];
 
-        // add area of CV and contribution from face with parent CV to diagonal
-        d[i] += cv_areas_[i] + a;
+        d[i] +=  a;
+        l[i]  = -a;
+        u[i]  = -a;
 
         // add contribution to the diagonal of parent
         d[p[i]] += a;
-
-        // note the symmetry
-        l[i] = -a;
-        u[i] = -a;
     }
 }
 
diff --git a/src/matrix.hpp b/src/matrix.hpp
index cd8fc30a9c18ebb90e846c0510cca7ac81855560..4fbd3346e274c1c712bc99d1bd4ea5c2d583f582 100644
--- a/src/matrix.hpp
+++ b/src/matrix.hpp
@@ -143,11 +143,13 @@ class matrix {
     void setup()
     {
         const auto n = size();
+        constexpr auto default_value
+            = std::numeric_limits<value_type>::quiet_NaN();
 
-        l_   = vector_type(n);
-        d_   = vector_type(n);
-        u_   = vector_type(n);
-        rhs_ = vector_type(n);
+        l_   = vector_type(n, default_value);
+        d_   = vector_type(n, default_value);
+        u_   = vector_type(n, default_value);
+        rhs_ = vector_type(n, default_value);
     }
 
     /// the parent indice that describe matrix structure
diff --git a/tests/test_run.cpp b/tests/test_run.cpp
index 6728cc20889792fe9e67883a4ab64777234bcd09..73e7929297f53c2972e52dd3116aecf0affa51ad 100644
--- a/tests/test_run.cpp
+++ b/tests/test_run.cpp
@@ -34,11 +34,25 @@ TEST(run, init)
     EXPECT_EQ(cell.soma()->mechanism("hh").get("gl").value, 0.0003);
     EXPECT_EQ(cell.soma()->mechanism("hh").get("el").value, -54.3);
 
-
-    cell.segment(1)->set_compartments(200);
+    cell.segment(1)->set_compartments(10);
 
     using fvm_cell = fvm::fvm_cell<double, int>;
     fvm_cell fvcell(cell);
+    EXPECT_EQ(fvcell.matrix().size(), 11);
+
+    fvcell.setup_matrx(0.01);
+
+    auto& J = fvcell.matrix();
+    //std::cout << "l" << J.l() << "\n";
+    //std::cout << "d" << J.d() << "\n";
+    //std::cout << "u" << J.u() << "\n";
+
+    J.rhs()(memory::all) = 1.;
+    J.rhs()[0] = 10.;
+
+    J.solve();
+
+    //std::cout << "x" << J.rhs() << "\n";
 }
 
 // test out the parameter infrastructure