diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp
index a2f6bd24eb..5b76a3edf3 100644
--- a/src/nrncvode/occvode.cpp
+++ b/src/nrncvode/occvode.cpp
@@ -1,3 +1,5 @@
+#include <Eigen/Eigen>
+
 #include <../../nrnconf.h>
 #include "hocdec.h"
 #include "cabcode.h"
@@ -17,10 +19,6 @@
 #include <cerrno>
 #include <numeric>
 
-
-#include "spmatrix.h"
-extern double* sp13mat;
-
 #if 1 || NRNMPI
 extern void (*nrnthread_v_transfer_)(NrnThread*);
 extern void (*nrnmpi_v_transfer_)();
@@ -349,7 +347,7 @@ void Cvode::daspk_init_eqn() {
     if (use_sparse13 == 0 || diam_changed != 0) {
         recalc_diam();
     }
-    zneq = spGetSize(_nt->_sp13mat, 0);
+    zneq = _nt->_sparseMat->cols();
     z.neq_v_ = z.nonvint_offset_ = zneq;
     // now add the membrane mechanism ode's to the count
     for (cml = z.cv_memb_list_; cml; cml = cml->next) {
diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp
index f07cfbed2b..96478b57db 100644
--- a/src/nrniv/matrixmap.cpp
+++ b/src/nrniv/matrixmap.cpp
@@ -2,7 +2,7 @@
 #include "matrixmap.h"
 #include <vector>
 
-#include "spmatrix.h"
+#include <Eigen/Eigen>
 
 MatrixMap::MatrixMap(Matrix& mat)
     : m_(mat) {}
@@ -19,7 +19,12 @@ void MatrixMap::mmfree() {
 
 void MatrixMap::add(double fac) {
     for (int i = 0; i < plen_; ++i) {
-        *ptree_[i] += fac * (*pm_[i]);
+        auto [j, k] = ptree_[i];
+        if (j != -1) {
+            // std::cerr << "modifying '" << nrn_threads->_sparseMat->coeff(j -1, k - 1) << "' (" <<
+            // j - 1 << ", " << k - 1 << ") with " << fac * (*pm_[i]) << std::endl;
+            nrn_threads->_sparseMat->coeffRef(j - 1, k - 1) += fac * (*pm_[i]);
+        }
     }
 }
 
@@ -27,11 +32,11 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) {
     NrnThread* _nt = nrn_threads;
     mmfree();
 
-    plen_ = 0;
     std::vector<int> nonzero_i, nonzero_j;
     m_.nonzeros(nonzero_i, nonzero_j);
-    pm_.resize(nonzero_i.size());
-    ptree_.resize(nonzero_i.size());
+    pm_.reserve(nonzero_i.size());
+    ptree_.reserve(nonzero_i.size());
+    plen_ = nonzero_i.size();
     for (int k = 0; k < nonzero_i.size(); k++) {
         const int i = nonzero_i[k];
         const int j = nonzero_j[k];
@@ -45,7 +50,7 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) {
             it = start + i - nnode;
         }
         int jt;
-        pm_[plen_] = m_.mep(i, j);
+        pm_.push_back(m_.mep(i, j));
         if (j < nnode) {
             jt = nodes[j]->eqn_index_ + layer[j];
             if (layer[j] > 0 && !nodes[j]->extnode) {
@@ -54,7 +59,13 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) {
         } else {
             jt = start + j - nnode;
         }
-        ptree_[plen_] = spGetElement(_nt->_sp13mat, it, jt);
-        ++plen_;
+        if (it == 0 || jt == 0) {
+            ptree_.emplace_back(-1, -1);
+        } else {
+            std::cerr << "someone touch me in the heart!: (" << it << ", " << jt << ")"
+                      << std::endl;
+            ptree_.emplace_back(it, jt);
+            _nt->_sparseMat->coeffRef(it - 1, jt - 1) = 0.;
+        }
     }
 }
diff --git a/src/nrniv/matrixmap.h b/src/nrniv/matrixmap.h
index f216dcc9c3..0a21e5ac0b 100644
--- a/src/nrniv/matrixmap.h
+++ b/src/nrniv/matrixmap.h
@@ -117,5 +117,5 @@ class MatrixMap {
     // the map
     int plen_ = 0;
     std::vector<double*> pm_{};
-    std::vector<double*> ptree_{};
+    std::vector<std::pair<int, int>> ptree_{};
 };
diff --git a/src/nrniv/nrndae.cpp b/src/nrniv/nrndae.cpp
index 4bbe44d8f6..2ab296ffab 100644
--- a/src/nrniv/nrndae.cpp
+++ b/src/nrniv/nrndae.cpp
@@ -8,7 +8,7 @@
 
 extern int secondorder;
 
-static NrnDAEPtrList nrndae_list;
+static std::list<NrnDAE*> nrndae_list;
 
 int nrndae_list_is_empty() {
     return nrndae_list.empty() ? 1 : 0;
@@ -221,7 +221,7 @@ void NrnDAE::dkmap(std::vector<neuron::container::data_handle<double>>& pv,
                                                                   y_.data() + i};
         pvdot[bmap_[i] - 1] =
             neuron::container::data_handle<double>{neuron::container::do_not_search,
-                                                   _nt->_sp13_rhs + bmap_[i]};
+                                                   _nt->_sparse_rhs + bmap_[i]};
     }
 }
 
@@ -231,7 +231,7 @@ void NrnDAE::update() {
     // note that the following is correct also for states that refer
     // to the internal potential of a segment. i.e rhs is v + vext[0]
     for (int i = 0; i < size_; ++i) {
-        y_[i] += _nt->_sp13_rhs[bmap_[i]];
+        y_[i] += _nt->_sparse_rhs[bmap_[i]];
     }
     // for (int i=0; i < size_; ++i) printf(" i=%d bmap_[i]=%d y_[i]=%g\n", i, bmap_[i],
     // y_->elem(i));
@@ -311,7 +311,7 @@ void NrnDAE::rhs() {
     v2y();
     f_(y_, yptmp_, size_);
     for (int i = 0; i < size_; ++i) {
-        _nt->_sp13_rhs[bmap_[i]] += yptmp_[i];
+        _nt->_sparse_rhs[bmap_[i]] += yptmp_[i];
     }
 }
 
diff --git a/src/nrnoc/fadvance.cpp b/src/nrnoc/fadvance.cpp
index a8a06f314e..305f0ec748 100644
--- a/src/nrnoc/fadvance.cpp
+++ b/src/nrnoc/fadvance.cpp
@@ -13,7 +13,7 @@
 #include "utils/profile/profiler_interface.h"
 #include "nonvintblock.h"
 #include "nrncvode.h"
-#include "spmatrix.h"
+#include <Eigen/Eigen>
 
 #include <vector>
 
@@ -672,11 +672,11 @@ void nrn_print_matrix(NrnThread* _nt) {
     Node* nd;
     if (use_sparse13) {
         if (ifarg(1) && chkarg(1, 0., 1.) == 0.) {
-            spPrint(_nt->_sp13mat, 1, 0, 1);
+            std::cout << &_nt->_sparseMat << std::endl;
         } else {
-            int i, n = spGetSize(_nt->_sp13mat, 0);
-            spPrint(_nt->_sp13mat, 1, 1, 1);
-            for (i = 1; i <= n; ++i) {
+            std::cout << &_nt->_sparseMat << std::endl;
+            int n = _nt->_sparseMat->cols();
+            for (int i = 1; i <= n; ++i) {
                 Printf("%d %g\n", i, _nt->actual_rhs(i));
             }
         }
diff --git a/src/nrnoc/multicore.cpp b/src/nrnoc/multicore.cpp
index b443e6c4b7..e347e2f8c6 100644
--- a/src/nrnoc/multicore.cpp
+++ b/src/nrnoc/multicore.cpp
@@ -48,6 +48,8 @@ the handling of v_structure_change as long as possible.
 #include <vector>
 #include <iostream>
 
+#include <Eigen/Eigen>
+
 #define CACHELINE_ALLOC(name, type, size) \
     name = (type*) nrn_cacheline_alloc((void**) &name, size * sizeof(type))
 #define CACHELINE_CALLOC(name, type, size) \
@@ -338,14 +340,14 @@ void nrn_threads_create(int n, bool parallel) {
                 for (j = 0; j < BEFORE_AFTER_SIZE; ++j) {
                     nt->tbl[j] = (NrnThreadBAList*) 0;
                 }
-                nt->_sp13_rhs = nullptr;
+                nt->_sparse_rhs = nullptr;
                 nt->_v_parent_index = 0;
                 nt->_v_node = 0;
                 nt->_v_parent = 0;
                 nt->_ecell_memb_list = 0;
                 nt->_ecell_child_cnt = 0;
                 nt->_ecell_children = NULL;
-                nt->_sp13mat = 0;
+                nt->_sparseMat = nullptr;
                 nt->_ctime = 0.0;
                 nt->_vcv = 0;
                 nt->_node_data_offset = 0;
@@ -439,9 +441,9 @@ void nrn_threads_free() {
             free(nt->_ecell_children);
             nt->_ecell_children = NULL;
         }
-        if (nt->_sp13mat) {
-            spDestroy(nt->_sp13mat);
-            nt->_sp13mat = 0;
+        if (nt->_sparseMat) {
+            delete nt->_sparseMat;
+            nt->_sparseMat = nullptr;
         }
         nt->end = 0;
         nt->ncell = 0;
diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h
index fb115b195c..5477a67cd7 100644
--- a/src/nrnoc/multicore.h
+++ b/src/nrnoc/multicore.h
@@ -27,6 +27,10 @@ actual_v, etc.
 */
 
 #include "membfunc.h"
+namespace Eigen {
+template <typename _Scalar, int _Options, typename _StorageIndex>
+class SparseMatrix;
+}
 
 #include <cstddef>
 
@@ -89,12 +93,12 @@ struct NrnThread {
     int* _v_parent_index;
     Node** _v_node;
     Node** _v_parent;
-    double* _sp13_rhs;           /* rhs matrix for sparse13 solver. updates to and from this vector
-                                    need to be transfered to actual_rhs */
-    char* _sp13mat;              /* handle to general sparse matrix */
-    Memb_list* _ecell_memb_list; /* normally nullptr */
-    Node** _ecell_children;      /* nodes with no extcell but parent has it */
-    void* _vcv;                  /* replaces old cvode_instance and nrn_cvode_ */
+    double* _sparse_rhs; /* rhs matrix for sparse13 solver. updates to and from this vector
+                          need to be transfered to actual_rhs */
+    Eigen::SparseMatrix<double, 1, int>* _sparseMat{}; /* handle to general sparse matrix */
+    Memb_list* _ecell_memb_list;                       /* normally nullptr */
+    Node** _ecell_children;                            /* nodes with no extcell but parent has it */
+    void* _vcv; /* replaces old cvode_instance and nrn_cvode_ */
 
 #if 1
     double _ctime; /* computation time in seconds (using nrnmpi_wtime) */
diff --git a/src/nrnoc/solve.cpp b/src/nrnoc/solve.cpp
index 38aa3ebc53..0a8029a531 100644
--- a/src/nrnoc/solve.cpp
+++ b/src/nrnoc/solve.cpp
@@ -57,7 +57,8 @@ node.v + extnode.v[0]
 #include "nrnmpiuse.h"
 #include "ocnotify.h"
 #include "section.h"
-#include "spmatrix.h"
+#include <Eigen/Eigen>
+#include <Eigen/LU>
 #include "treeset.h"
 
 #include <stdio.h>
@@ -364,22 +365,13 @@ void nrn_solve(NrnThread* _nt) {
     }
 #else
     if (use_sparse13) {
-        int e;
         nrn_thread_error("solve use_sparse13");
         update_sp13_mat_based_on_actual_d(_nt);
-        e = spFactor(_nt->_sp13mat);
-        if (e != spOKAY) {
-            switch (e) {
-            case spZERO_DIAG:
-                hoc_execerror("spFactor error:", "Zero Diagonal");
-            case spNO_MEMORY:
-                hoc_execerror("spFactor error:", "No Memory");
-            case spSINGULAR:
-                hoc_execerror("spFactor error:", "Singular");
-            }
-        }
         update_sp13_rhs_based_on_actual_rhs(_nt);
-        spSolve(_nt->_sp13mat, _nt->_sp13_rhs, _nt->_sp13_rhs);
+        _nt->_sparseMat->makeCompressed();
+        Eigen::SparseLU<Eigen::SparseMatrix<double, Eigen::RowMajor>> lu(*_nt->_sparseMat);
+        Eigen::Map<Eigen::VectorXd> rhs(_nt->_sparse_rhs + 1, _nt->_sparseMat->cols());
+        rhs = lu.solve(rhs);
         update_actual_d_based_on_sp13_mat(_nt);
         update_actual_rhs_based_on_sp13_rhs(_nt);
     } else {
diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp
index 2a1c0abf97..f3637addb0 100644
--- a/src/nrnoc/treeset.cpp
+++ b/src/nrnoc/treeset.cpp
@@ -18,7 +18,7 @@
 #include "ocnotify.h"
 #include "partrans.h"
 #include "section.h"
-#include "spmatrix.h"
+#include <Eigen/Eigen>
 #include "utils/profile/profiler_interface.h"
 #include "multicore.h"
 
@@ -32,8 +32,6 @@
 
 #include <fmt/format.h>
 
-extern spREAL* spGetElement(char*, int, int);
-
 int nrn_shape_changed_; /* for notifying Shape class in nrniv */
 double* nrn_mech_wtime_;
 
@@ -331,20 +329,20 @@ vm += dvi-dvx
 */
 
 /*
- * Update actual_rhs based on _sp13_rhs used for sparse13 solver
+ * Update actual_rhs based on _sparse_rhs used for sparse13 solver
  */
 void update_actual_rhs_based_on_sp13_rhs(NrnThread* nt) {
     for (int i = 0; i < nt->end; i++) {
-        nt->actual_rhs(i) = nt->_sp13_rhs[nt->_v_node[i]->eqn_index_];
+        nt->actual_rhs(i) = nt->_sparse_rhs[nt->_v_node[i]->eqn_index_];
     }
 }
 
 /*
- * Update _sp13_rhs used for sparse13 solver based on changes on actual_rhs
+ * Update _sparse_rhs used for sparse13 solver based on changes on actual_rhs
  */
 void update_sp13_rhs_based_on_actual_rhs(NrnThread* nt) {
     for (int i = 0; i < nt->end; i++) {
-        nt->_sp13_rhs[nt->_v_node[i]->eqn_index_] = nt->actual_rhs(i);
+        nt->_sparse_rhs[nt->_v_node[i]->eqn_index_] = nt->actual_rhs(i);
     }
 }
 
@@ -392,13 +390,12 @@ void nrn_rhs(neuron::model_sorted_token const& cache_token, NrnThread& nt) {
     }
     auto* const vec_rhs = nt.node_rhs_storage();
     if (use_sparse13) {
-        int i, neqn;
         nrn_thread_error("nrn_rhs use_sparse13");
-        neqn = spGetSize(_nt->_sp13mat, 0);
-        for (i = 1; i <= neqn; ++i) {
-            _nt->_sp13_rhs[i] = 0.;
+        int neqn = _nt->_sparseMat->cols();
+        for (int i = 1; i <= neqn; ++i) {
+            _nt->_sparse_rhs[i] = 0.;
         }
-        for (i = i1; i < i3; ++i) {
+        for (int i = i1; i < i3; ++i) {
             NODERHS(_nt->_v_node[i]) = 0.;
         }
     } else {
@@ -499,8 +496,12 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) {
     }
 
     if (use_sparse13) {
-        // Zero the sparse13 matrix
-        spClear(_nt->_sp13mat);
+        Eigen::SparseMatrix<double, Eigen::RowMajor>& m_ = *_nt->_sparseMat;
+        for (int k = 0; k < m_.outerSize(); ++k) {
+            for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(m_, k); it; ++it) {
+                it.valueRef() = 0.;
+            }
+        }
     }
 
     // Make sure the SoA node diagonals are also zeroed (is this needed?)
@@ -1785,12 +1786,12 @@ void v_setup_vectors(void) {
 void nrn_matrix_node_free() {
     NrnThread* nt;
     FOR_THREADS(nt) {
-        if (nt->_sp13_rhs) {
-            free(std::exchange(nt->_sp13_rhs, nullptr));
+        if (nt->_sparse_rhs) {
+            free(std::exchange(nt->_sparse_rhs, nullptr));
         }
-        if (nt->_sp13mat) {
-            spDestroy(nt->_sp13mat);
-            nt->_sp13mat = (char*) 0;
+        if (nt->_sparseMat) {
+            delete nt->_sparseMat;
+            nt->_sparseMat = nullptr;
         }
     }
     diam_changed = 1;
@@ -1853,14 +1854,6 @@ int nrn_method_consistent(void) {
 }
 
 
-/*
-sparse13 uses the mathematical index scheme in which the rows and columns
-range from 1...size instead of 0...size-1. Thus the calls to
-spGetElement(i,j) have args that are 1 greater than a normal c style.
-Also the actual_rhs_ uses this style, 1-neqn, when sparse13 is activated
-and therefore is passed to spSolve as actual_rhs intead of actual_rhs-1.
-*/
-
 static void nrn_matrix_node_alloc(void) {
     int i, b;
     Node* nd;
@@ -1868,16 +1861,16 @@ static void nrn_matrix_node_alloc(void) {
 
     nrn_method_consistent();
     nt = nrn_threads;
-    /*printf("use_sparse13=%d sp13mat=%lx rhs=%lx\n", use_sparse13, (long)nt->_sp13mat,
+    /*printf("use_sparse13=%d sp13mat=%lx rhs=%lx\n", use_sparse13, (long)nt->_sparseMat,
      * (long)nt->_actual_rhs);*/
     if (use_sparse13) {
-        if (nt->_sp13mat) {
+        if (nt->_sparseMat) {
             return;
         } else {
             nrn_matrix_node_free();
         }
     } else {
-        if (nt->_sp13mat) {
+        if (nt->_sparseMat) {
             v_structure_change = 1;
             v_setup_vectors();
             return;
@@ -1887,7 +1880,7 @@ static void nrn_matrix_node_alloc(void) {
     }
     ++nrn_matrix_cnt_;
     if (use_sparse13) {
-        int in, err, extn, neqn, j;
+        int in, extn, neqn, j;
         nt = nrn_threads;
         neqn = nt->end + nrndae_extra_eqn_count();
         extn = 0;
@@ -1896,53 +1889,94 @@ static void nrn_matrix_node_alloc(void) {
         }
         /*printf(" %d extracellular nodes\n", extn);*/
         neqn += extn;
-        nt->_sp13_rhs = (double*) ecalloc(neqn + 1, sizeof(double));
-        nt->_sp13mat = spCreate(neqn, 0, &err);
-        if (err != spOKAY) {
-            hoc_execerror("Couldn't create sparse matrix", (char*) 0);
-        }
+        nt->_sparse_rhs = (double*) ecalloc(neqn + 1, sizeof(double));
+        nt->_sparseMat = new Eigen::SparseMatrix<double, Eigen::RowMajor>(neqn, neqn);
         for (in = 0, i = 1; in < nt->end; ++in, ++i) {
             nt->_v_node[in]->eqn_index_ = i;
             if (nt->_v_node[in]->extnode) {
                 i += nlayer;
             }
         }
-        for (in = 0; in < nt->end; ++in) {
-            int ie, k;
-            Node *nd, *pnd;
-            Extnode* nde;
-            nd = nt->_v_node[in];
-            nde = nd->extnode;
-            pnd = nt->_v_parent[in];
-            i = nd->eqn_index_;
-            nt->_sp13_rhs[i] = nt->actual_rhs(in);
-            nd->_d_matelm = spGetElement(nt->_sp13mat, i, i);
-            if (nde) {
-                for (ie = 0; ie < nlayer; ++ie) {
-                    k = i + ie + 1;
-                    nde->_d[ie] = spGetElement(nt->_sp13mat, k, k);
-                    nde->_rhs[ie] = nt->_sp13_rhs + k;
-                    nde->_x21[ie] = spGetElement(nt->_sp13mat, k, k - 1);
-                    nde->_x12[ie] = spGetElement(nt->_sp13mat, k - 1, k);
+        {
+            // sparse13 was a linked list of single values and so we were able to have
+            // stable pointers on internal values.
+            // Eigen::SparseMatrix is more like a vector, that is reallocated when needed
+            // and so pointers will be invalidate.
+            // More over, when solving with SparseLU, the matrix is compressed and the
+            // internal state is fully changed.
+            // setFromTriplets is an efficient way of having an already compressed matrix.
+            // We can safely get stable pointers as long as we don't insert new value with
+            // coeff() and coeffRef().
+            std::vector<Eigen::Triplet<double>> triplets{};
+            for (in = 0; in < nt->end; ++in) {
+                const Node* nd = nt->_v_node[in];
+                const Extnode* nde = nd->extnode;
+                const Node* pnd = nt->_v_parent[in];
+                int i = nd->eqn_index_;
+                triplets.emplace_back(i - 1, i - 1, 0.);
+                if (nde) {
+                    for (int ie = 0; ie < nlayer; ++ie) {
+                        int k = i + ie + 1;
+                        triplets.emplace_back(k - 1, k - 1, 0.);
+                        triplets.emplace_back(k - 1, k - 2, 0.);
+                        triplets.emplace_back(k - 2, k - 1, 0.);
+                    }
+                }
+                if (pnd) {
+                    int j = pnd->eqn_index_;
+                    triplets.emplace_back(j - 1, i - 1, 0.);
+                    triplets.emplace_back(i - 1, j - 1, 0.);
+                    if (nde && pnd->extnode)
+                        for (int ie = 0; ie < nlayer; ++ie) {
+                            int kp = j + ie + 1;
+                            int k = i + ie + 1;
+                            triplets.emplace_back(kp - 1, k - 1, 0.);
+                            triplets.emplace_back(k - 1, kp - 1, 0.);
+                        }
                 }
             }
-            if (pnd) {
-                j = pnd->eqn_index_;
-                nd->_a_matelm = spGetElement(nt->_sp13mat, j, i);
-                nd->_b_matelm = spGetElement(nt->_sp13mat, i, j);
-                if (nde && pnd->extnode)
-                    for (ie = 0; ie < nlayer; ++ie) {
-                        int kp = j + ie + 1;
-                        k = i + ie + 1;
-                        nde->_a_matelm[ie] = spGetElement(nt->_sp13mat, kp, k);
-                        nde->_b_matelm[ie] = spGetElement(nt->_sp13mat, k, kp);
+            nt->_sparseMat->setFromTriplets(triplets.begin(), triplets.end());
+        }
+        // Add some new nodes
+        nrndae_alloc();
+        nt->_sparseMat->makeCompressed();
+        {
+            // Now that we have our stable matrix, let's extract pointers.
+            for (in = 0; in < nt->end; ++in) {
+                Node *nd, *pnd;
+                Extnode* nde;
+                nd = nt->_v_node[in];
+                nde = nd->extnode;
+                pnd = nt->_v_parent[in];
+                i = nd->eqn_index_;
+                nt->_sparse_rhs[i] = nt->actual_rhs(in);
+                nd->_d_matelm = &nt->_sparseMat->coeffRef(i - 1, i - 1);
+                if (nde) {
+                    for (int ie = 0; ie < nlayer; ++ie) {
+                        int k = i + ie + 1;
+                        nde->_d[ie] = &nt->_sparseMat->coeffRef(k - 1, k - 1);
+                        nde->_rhs[ie] = nt->_sparse_rhs + k;
+                        nde->_x21[ie] = &nt->_sparseMat->coeffRef(k - 1, k - 2);
+                        nde->_x12[ie] = &nt->_sparseMat->coeffRef(k - 2, k - 1);
                     }
-            } else { /* not needed if index starts at 1 */
-                nd->_a_matelm = nullptr;
-                nd->_b_matelm = nullptr;
+                }
+                if (pnd) {
+                    j = pnd->eqn_index_;
+                    nd->_a_matelm = &nt->_sparseMat->coeffRef(j - 1, i - 1);
+                    nd->_b_matelm = &nt->_sparseMat->coeffRef(i - 1, j - 1);
+                    if (nde && pnd->extnode)
+                        for (int ie = 0; ie < nlayer; ++ie) {
+                            int kp = j + ie + 1;
+                            int k = i + ie + 1;
+                            nde->_a_matelm[ie] = &nt->_sparseMat->coeffRef(kp - 1, k - 1);
+                            nde->_b_matelm[ie] = &nt->_sparseMat->coeffRef(k - 1, kp - 1);
+                        }
+                } else { /* not needed if index starts at 1 */
+                    nd->_a_matelm = nullptr;
+                    nd->_b_matelm = nullptr;
+                }
             }
         }
-        nrndae_alloc();
     } else {
         FOR_THREADS(nt) {
             assert(nrndae_extra_eqn_count() == 0);