diff --git a/src/nocmodl/noccout.cpp b/src/nocmodl/noccout.cpp
index 31c83e2e73..9a7952b06d 100644
--- a/src/nocmodl/noccout.cpp
+++ b/src/nocmodl/noccout.cpp
@@ -260,7 +260,8 @@ void c_out() {
             P("  }\n");
             P("#if EXTRACELLULAR\n");
             P(" if (auto* const _extnode = _nrn_mechanism_access_extnode(_nd); _extnode) {\n");
-            P("   *_extnode->_d[0] += _g;\n");
+            P("   int ind = _nrn_mechanism_access_index(_nd);\n");
+            P("   _nrn_mechanism_get_matrix_elem(_nt, ind, ind) += _g;\n");
             P(" }\n");
             P("#endif\n");
         } else {
@@ -679,7 +680,8 @@ void c_out_vectorize() {
             P("  }\n");
             P("#if EXTRACELLULAR\n");
             P(" if (auto* const _extnode = _nrn_mechanism_access_extnode(_nd); _extnode) {\n");
-            P("   *_extnode->_d[0] += _g;\n");
+            P("   int ind = _nrn_mechanism_access_index(_nd);\n");
+            P("   _nrn_mechanism_get_matrix_elem(_nt, ind, ind) += _g;\n");
             P(" }\n");
             P("#endif\n");
         } else {
diff --git a/src/nrncvode/nrndaspk.cpp b/src/nrncvode/nrndaspk.cpp
index 5434cdc031..255ca35733 100644
--- a/src/nrncvode/nrndaspk.cpp
+++ b/src/nrncvode/nrndaspk.cpp
@@ -6,7 +6,6 @@
 
 #include <stdio.h>
 #include <math.h>
-#include "spmatrix.h"
 #include "nrnoc2iv.h"
 #include "cvodeobj.h"
 #include "nrndaspk.h"
diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp
index a2f6bd24eb..e81513a61c 100644
--- a/src/nrncvode/occvode.cpp
+++ b/src/nrncvode/occvode.cpp
@@ -18,7 +18,7 @@
 #include <numeric>
 
 
-#include "spmatrix.h"
+#include "ocmatrix.h"
 extern double* sp13mat;
 
 #if 1 || NRNMPI
@@ -349,7 +349,7 @@ void Cvode::daspk_init_eqn() {
     if (use_sparse13 == 0 || diam_changed != 0) {
         recalc_diam();
     }
-    zneq = spGetSize(_nt->_sp13mat, 0);
+    zneq = _nt->_sp13mat->ncol();
     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 cc584e0942..69d3a85249 100644
--- a/src/nrniv/matrixmap.cpp
+++ b/src/nrniv/matrixmap.cpp
@@ -2,8 +2,6 @@
 #include "matrixmap.h"
 #include <vector>
 
-#include "spmatrix.h"
-
 MatrixMap::MatrixMap(Matrix& mat)
     : m_(mat) {}
 
@@ -18,9 +16,11 @@ void MatrixMap::mmfree() {
 }
 
 void MatrixMap::add(double fac) {
+    NrnThread* _nt = nrn_threads;
     for (int i = 0; i < pm_.size(); ++i) {
         auto [it, jt] = pm_[i];
-        *ptree_[i] += fac * m_(it, jt);
+        auto [im, jm] = ptree_[i];
+        _nt->_sp13mat->coeff(im, jm) += fac * m_(it, jt);
     }
 }
 
@@ -38,7 +38,6 @@ int MatrixMap::compute_index(int i, int start, int nnode, Node** nodes, int* lay
 }
 
 void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) {
-    NrnThread* _nt = nrn_threads;
     mmfree();
 
     std::vector<std::pair<int, int>> nzs = m_.nonzeros();
@@ -49,7 +48,7 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) {
         int jt = compute_index(j, start, nnode, nodes, layer);
         if (it != 0 && jt != 0) {
             pm_.emplace_back(i, j);
-            ptree_.emplace_back(spGetElement(_nt->_sp13mat, it, jt));
+            ptree_.emplace_back(it - 1, jt - 1);
         }
     }
 }
diff --git a/src/nrniv/matrixmap.h b/src/nrniv/matrixmap.h
index c234e8369c..beffa15d9a 100644
--- a/src/nrniv/matrixmap.h
+++ b/src/nrniv/matrixmap.h
@@ -116,7 +116,7 @@ class MatrixMap {
 
     // the map
     std::vector<std::pair<int, int>> pm_{};
-    std::vector<double*> ptree_{};
+    std::vector<std::pair<int, int>> ptree_{};
 
     int compute_index(int, int, int, Node**, int*) const;
 };
diff --git a/src/nrnoc/extcelln.cpp b/src/nrnoc/extcelln.cpp
index b6facf2b76..552e9b58cb 100644
--- a/src/nrnoc/extcelln.cpp
+++ b/src/nrnoc/extcelln.cpp
@@ -7,6 +7,7 @@
 #include "nrniv_mf.h"
 #include "hocassrt.h"
 #include "parse.hpp"
+#include "ocmatrix.h"
 
 
 extern int nrn_use_daspk_;
@@ -213,16 +214,11 @@ static void extcell_init(neuron::model_sorted_token const&,
 void extnode_free_elements(Extnode* nde) {
     if (nde->v) {
         free(nde->v);  /* along with _a and _b */
-        free(nde->_d); /* along with _rhs, _a_matelm, _b_matelm, _x12, and _x21 */
-        nde->v = NULL;
-        nde->_a = NULL;
-        nde->_b = NULL;
-        nde->_d = NULL;
-        nde->_rhs = NULL;
-        nde->_a_matelm = NULL;
-        nde->_b_matelm = NULL;
-        nde->_x12 = NULL;
-        nde->_x21 = NULL;
+        free(nde->_rhs);
+        nde->v = nullptr;
+        nde->_a = nullptr;
+        nde->_b = nullptr;
+        nde->_rhs = nullptr;
     }
 }
 
@@ -294,12 +290,7 @@ static void extnode_alloc_elements(Extnode* nde) {
         nde->_a = nde->v + nlayer;
         nde->_b = nde->_a + nlayer;
 
-        nde->_d = (double**) ecalloc(nlayer * 6, sizeof(double*));
-        nde->_rhs = nde->_d + nlayer;
-        nde->_a_matelm = nde->_rhs + nlayer;
-        nde->_b_matelm = nde->_a_matelm + nlayer;
-        nde->_x12 = nde->_b_matelm + nlayer;
-        nde->_x21 = nde->_x12 + nlayer;
+        nde->_rhs = (double**) ecalloc(nlayer * 1, sizeof(double*));
     }
 }
 
@@ -426,7 +417,6 @@ void nrn_rhs_ext(NrnThread* _nt) {
 void nrn_setup_ext(NrnThread* _nt) {
     int i, j, cnt;
     Node *nd, *pnd, **ndlist;
-    double* pd;
     double d, cfac, mfac;
     Extnode *nde, *pnde;
     Memb_list* ml = _nt->_ecell_memb_list;
@@ -441,47 +431,56 @@ void nrn_setup_ext(NrnThread* _nt) {
     /* d contains all the membrane conductances (and capacitance) */
     /* i.e. (cm/dt + di/dvm - dis/dvi)*[dvi] and
         (dis/dvi)*[dvx] */
+    // This loop handle conductances between the node and the first layer
     for (i = 0; i < cnt; ++i) {
+        OcMatrix& m = *_nt->_sp13mat;
         nd = ndlist[i];
+        int index = nd->eqn_index_;
         nde = nd->extnode;
+        int ext_index = nde->eqn_index_;
         d = NODED(nd);
         /* nde->_d only has -ELECTRODE_CURRENT contribution */
-        d = (*nde->_d[0] += NODED(nd));
+        m(ext_index + 0, ext_index + 0) += NODED(nd);
+        d = m.getval(ext_index + 0, ext_index + 0);
         /* now d is only the membrane current contribution */
         /* i.e. d =  cm/dt + di/dvm */
-        *nde->_x12[0] -= d;
-        *nde->_x21[0] -= d;
+        m(index - 1, ext_index + 0) -= d;
+        m(ext_index + 0, index - 1) -= d;
 #if I_MEMBRANE
         ml->data(i, sav_g_index) = d;
 #endif
     }
     /* series resistance, capacitance, and axial terms. */
+    // This look takes care of the conductances between layers
     for (i = 0; i < cnt; ++i) {
+        OcMatrix& m = *_nt->_sp13mat;
         nd = ndlist[i];
         nde = nd->extnode;
+        int ext_index = nde->eqn_index_;
         pnd = _nt->_v_parent[nd->v_node_index];
         if (pnd) {
             /* series resistance and capacitance to ground */
             j = 0;
             for (;;) { /* between j and j+1 layer */
                 mfac = (*nde->param[xg_index_ext(j)] + *nde->param[xc_index_ext(j)] * cfac);
-                *nde->_d[j] += mfac;
+                m(ext_index + j, ext_index + j) += mfac;
                 ++j;
                 if (j == nrn_nlayer_extracellular) {
                     break;
                 }
-                *nde->_d[j] += mfac;
-                *nde->_x12[j] -= mfac;
-                *nde->_x21[j] -= mfac;
+                m(ext_index + j, ext_index + j) += mfac;
+                m(ext_index + j - 1, ext_index + j) -= mfac;
+                m(ext_index + j, ext_index + j - 1) -= mfac;
             }
             pnde = pnd->extnode;
             /* axial connections */
             if (pnde) { /* parent sec may not be extracellular */
+                int parent_ext_index = pnde->eqn_index_;
                 for (j = 0; j < nrn_nlayer_extracellular; ++j) {
-                    *nde->_d[j] -= nde->_b[j];
-                    *pnde->_d[j] -= nde->_a[j];
-                    *nde->_a_matelm[j] += nde->_a[j];
-                    *nde->_b_matelm[j] += nde->_b[j];
+                    m(ext_index + j, ext_index + j) -= nde->_b[j];
+                    m(parent_ext_index + j, parent_ext_index + j) -= nde->_a[j];
+                    m(parent_ext_index + j, ext_index + j) += nde->_a[j];
+                    m(ext_index + j, parent_ext_index + j) += nde->_b[j];
                 }
             }
         }
diff --git a/src/nrnoc/fadvance.cpp b/src/nrnoc/fadvance.cpp
index a8a06f314e..82c9e901c9 100644
--- a/src/nrnoc/fadvance.cpp
+++ b/src/nrnoc/fadvance.cpp
@@ -13,7 +13,6 @@
 #include "utils/profile/profiler_interface.h"
 #include "nonvintblock.h"
 #include "nrncvode.h"
-#include "spmatrix.h"
 
 #include <vector>
 
@@ -672,13 +671,13 @@ 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);
+            // spPrint(_nt->_sp13mat, 1, 0, 1);
         } else {
-            int i, n = spGetSize(_nt->_sp13mat, 0);
-            spPrint(_nt->_sp13mat, 1, 1, 1);
-            for (i = 1; i <= n; ++i) {
-                Printf("%d %g\n", i, _nt->actual_rhs(i));
-            }
+            // int i, n = spGetSize(_nt->_sp13mat, 0);
+            // spPrint(_nt->_sp13mat, 1, 1, 1);
+            // for (i = 1; i <= n; ++i) {
+            //     Printf("%d %g\n", i, _nt->actual_rhs(i));
+            // }
         }
     } else if (_nt) {
         for (inode = 0; inode < _nt->end; ++inode) {
diff --git a/src/nrnoc/membfunc.cpp b/src/nrnoc/membfunc.cpp
index e858536147..7d48a0ee0b 100644
--- a/src/nrnoc/membfunc.cpp
+++ b/src/nrnoc/membfunc.cpp
@@ -2,6 +2,7 @@
 
 #include "multicore.h"
 #include "section.h"
+#include "ocmatrix.h"
 
 #include <cassert>
 
@@ -49,6 +50,13 @@ double& _nrn_mechanism_access_rhs(Node* node) {
 double& _nrn_mechanism_access_voltage(Node* node) {
     return node->v();
 }
+int _nrn_mechanism_access_index(const Node* node) {
+    return node->eqn_index_;
+}
+double& _nrn_mechanism_get_matrix_elem(NrnThread* nt, int i, int j) {
+    OcMatrix& m = *nt->_sp13mat;
+    return m.coeff(i, j);
+}
 neuron::container::data_handle<double> _nrn_mechanism_get_area_handle(Node* node) {
     if (node) {
         return node->area_handle();
diff --git a/src/nrnoc/membfunc.h b/src/nrnoc/membfunc.h
index 18647c19b7..5d7d79b8be 100644
--- a/src/nrnoc/membfunc.h
+++ b/src/nrnoc/membfunc.h
@@ -306,6 +306,8 @@ namespace _get {
 [[nodiscard]] double& _nrn_mechanism_access_param(Prop*, int field, int array_index = 0);
 [[nodiscard]] double& _nrn_mechanism_access_rhs(Node*);
 [[nodiscard]] double& _nrn_mechanism_access_voltage(Node*);
+[[nodiscard]] int _nrn_mechanism_access_index(const Node*);
+[[nodiscard]] double& _nrn_mechanism_get_matrix_elem(NrnThread* nt, int, int );
 [[nodiscard]] neuron::container::data_handle<double> _nrn_mechanism_get_area_handle(Node*);
 [[nodiscard]] Section* _nrn_mechanism_get_child(Section*);
 [[nodiscard]] int _nrn_mechanism_get_nnode(Section*);
diff --git a/src/nrnoc/multicore.cpp b/src/nrnoc/multicore.cpp
index b443e6c4b7..9c01395993 100644
--- a/src/nrnoc/multicore.cpp
+++ b/src/nrnoc/multicore.cpp
@@ -38,6 +38,7 @@ the handling of v_structure_change as long as possible.
 */
 
 #include "nmodlmutex.h"
+#include "ocmatrix.h"
 
 #include <cstdint>
 #include <condition_variable>
@@ -345,7 +346,7 @@ void nrn_threads_create(int n, bool parallel) {
                 nt->_ecell_memb_list = 0;
                 nt->_ecell_child_cnt = 0;
                 nt->_ecell_children = NULL;
-                nt->_sp13mat = 0;
+                nt->_sp13mat = nullptr;
                 nt->_ctime = 0.0;
                 nt->_vcv = 0;
                 nt->_node_data_offset = 0;
@@ -440,8 +441,8 @@ void nrn_threads_free() {
             nt->_ecell_children = NULL;
         }
         if (nt->_sp13mat) {
-            spDestroy(nt->_sp13mat);
-            nt->_sp13mat = 0;
+            delete(nt->_sp13mat);
+            nt->_sp13mat = nullptr;
         }
         nt->end = 0;
         nt->ncell = 0;
diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h
index fb115b195c..0511fbb33b 100644
--- a/src/nrnoc/multicore.h
+++ b/src/nrnoc/multicore.h
@@ -29,6 +29,7 @@ actual_v, etc.
 #include "membfunc.h"
 
 #include <cstddef>
+class OcMatrix;
 
 typedef struct NrnThreadMembList { /* patterned after CvMembList in cvodeobj.h */
     struct NrnThreadMembList* next;
@@ -91,7 +92,7 @@ struct NrnThread {
     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 */
+    OcMatrix* _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_ */
diff --git a/src/nrnoc/section.h b/src/nrnoc/section.h
index 1e26dc8824..760b2b8764 100644
--- a/src/nrnoc/section.h
+++ b/src/nrnoc/section.h
@@ -181,9 +181,6 @@ struct Node {
         return _node_handle.non_owning_handle();
     }
     double _rinv{}; /* conductance uS from node to parent */
-    double* _a_matelm;
-    double* _b_matelm;
-    double* _d_matelm;
     int eqn_index_;                 /* sparse13 matrix row/col index */
                                     /* if no extnodes then = v_node_index +1*/
                                     /* each extnode adds nlayer more equations after this */
diff --git a/src/nrnoc/section_fwd.hpp b/src/nrnoc/section_fwd.hpp
index fe1a955f6c..ee86737e13 100644
--- a/src/nrnoc/section_fwd.hpp
+++ b/src/nrnoc/section_fwd.hpp
@@ -40,12 +40,9 @@ struct Extnode {
     double* v; /* v external. */
     double* _a;
     double* _b;
-    double** _d;
-    double** _rhs; /* d, rhs, a, and b are analogous to those in node */
-    double** _a_matelm;
-    double** _b_matelm;
-    double** _x12; /* effect of v[layer] on eqn layer-1 (or internal)*/
-    double** _x21; /* effect of v[layer-1 or internal] on eqn layer*/
+    double** _rhs; /* rhs, a, and b are analogous to those in node */
+
+    int eqn_index_;
 };
 #endif
 
diff --git a/src/nrnoc/solve.cpp b/src/nrnoc/solve.cpp
index 38aa3ebc53..8c4b0286af 100644
--- a/src/nrnoc/solve.cpp
+++ b/src/nrnoc/solve.cpp
@@ -57,8 +57,9 @@ node.v + extnode.v[0]
 #include "nrnmpiuse.h"
 #include "ocnotify.h"
 #include "section.h"
-#include "spmatrix.h"
+#include "ocmatrix.h"
 #include "treeset.h"
+#include "ivocvect.h"
 
 #include <stdio.h>
 #include <stdlib.h>
@@ -367,19 +368,17 @@ void nrn_solve(NrnThread* _nt) {
         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);
+        int size = _nt->_sp13mat->ncol();
+        std::vector<double> vector_in(_nt->_sp13_rhs + 1, _nt->_sp13_rhs + size + 1);
+        IvocVect in; in.vec() = vector_in;
+        IvocVect out(size);
+        _nt->_sp13mat->solv(&in, &out, false);
+        for (int i = 0; i < size; ++i) {
+            _nt->_sp13_rhs[i + 1] = out[i];
+
+        }
         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..cfb3f3c2d0 100644
--- a/src/nrnoc/treeset.cpp
+++ b/src/nrnoc/treeset.cpp
@@ -18,9 +18,10 @@
 #include "ocnotify.h"
 #include "partrans.h"
 #include "section.h"
-#include "spmatrix.h"
+#include "ocmatrix.h"
 #include "utils/profile/profiler_interface.h"
 #include "multicore.h"
+#include "ocmatrix.h"
 
 #include <cstdio>
 #include <cstdlib>
@@ -32,8 +33,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_;
 
@@ -353,7 +352,9 @@ void update_sp13_rhs_based_on_actual_rhs(NrnThread* nt) {
  */
 void update_actual_d_based_on_sp13_mat(NrnThread* nt) {
     for (int i = 0; i < nt->end; ++i) {
-        nt->actual_d(i) = *nt->_v_node[i]->_d_matelm;
+        const int index = nt->_v_node[i]->eqn_index_;
+        OcMatrix& m = *nt->_sp13mat;
+        nt->actual_d(i) = m.getval(index - 1, index - 1);
     }
 }
 
@@ -362,7 +363,9 @@ void update_actual_d_based_on_sp13_mat(NrnThread* nt) {
  */
 void update_sp13_mat_based_on_actual_d(NrnThread* nt) {
     for (int i = 0; i < nt->end; ++i) {
-        *nt->_v_node[i]->_d_matelm = nt->actual_d(i);
+        const int index = nt->_v_node[i]->eqn_index_;
+        OcMatrix& m = *nt->_sp13mat;
+        m(index - 1, index - 1) = nt->actual_d(i);
     }
 }
 
@@ -394,7 +397,7 @@ void nrn_rhs(neuron::model_sorted_token const& cache_token, NrnThread& nt) {
     if (use_sparse13) {
         int i, neqn;
         nrn_thread_error("nrn_rhs use_sparse13");
-        neqn = spGetSize(_nt->_sp13mat, 0);
+        neqn = _nt->_sp13mat->nrow();
         for (i = 1; i <= neqn; ++i) {
             _nt->_sp13_rhs[i] = 0.;
         }
@@ -500,7 +503,7 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) {
 
     if (use_sparse13) {
         // Zero the sparse13 matrix
-        spClear(_nt->_sp13mat);
+        _nt->_sp13mat->zero();
     }
 
     // Make sure the SoA node diagonals are also zeroed (is this needed?)
@@ -550,7 +553,6 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) {
         }
     }
 #if EXTRACELLULAR
-    /* nde->_d[0] contains the -ELECTRODE_CURRENT contribution to nd->_d */
     nrn_setup_ext(_nt);
 #endif
     if (use_sparse13) {
@@ -571,16 +573,21 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) {
         update_sp13_mat_based_on_actual_d(_nt);  // just because of activclamp_lhs
         for (i = i2; i < i3; ++i) {              // note i2
             Node* nd = _nt->_v_node[i];
+            OcMatrix& m = *_nt->_sp13mat;
             auto const parent_i = _nt->_v_parent_index[i];
             auto* const parent_nd = _nt->_v_node[parent_i];
             auto const nd_a = NODEA(nd);
             auto const nd_b = NODEB(nd);
             // Update entries in sp13_mat
-            *nd->_a_matelm += nd_a;
-            *nd->_b_matelm += nd_b; /* b may have value from lincir */
-            *nd->_d_matelm -= nd_b;
-            // used to update NODED (sparse13 matrix) using NODEA and NODEB ("SoA")
-            *parent_nd->_d_matelm -= nd_a;
+            {
+                const int index = nd->eqn_index_;
+                const int parent_index = parent_nd->eqn_index_;
+                m(parent_index - 1, index - 1) += nd_a;
+                m(index - 1, parent_index - 1) += nd_b; /* b may have value from lincir */
+                m(index - 1, index - 1) -= nd_b;
+                // used to update NODED (sparse13 matrix) using NODEA and NODEB ("SoA")
+                m(parent_index - 1, parent_index - 1) -= nd_a;
+            }
             // Also update the Node's d value in the SoA storage (is this needed?)
             vec_d[i] -= nd_b;
             vec_d[parent_i] -= nd_a;
@@ -1789,8 +1796,8 @@ void nrn_matrix_node_free() {
             free(std::exchange(nt->_sp13_rhs, nullptr));
         }
         if (nt->_sp13mat) {
-            spDestroy(nt->_sp13mat);
-            nt->_sp13mat = (char*) 0;
+            delete nt->_sp13mat;
+            nt->_sp13mat = nullptr;
         }
     }
     diam_changed = 1;
@@ -1887,7 +1894,7 @@ static void nrn_matrix_node_alloc(void) {
     }
     ++nrn_matrix_cnt_;
     if (use_sparse13) {
-        int in, err, extn, neqn, j;
+        int in, extn, neqn;
         nt = nrn_threads;
         neqn = nt->end + nrndae_extra_eqn_count();
         extn = 0;
@@ -1897,50 +1904,28 @@ 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->_sp13mat = OcMatrix::instance(neqn, neqn, OcMatrix::MSPARSE);
         for (in = 0, i = 1; in < nt->end; ++in, ++i) {
-            nt->_v_node[in]->eqn_index_ = i;
+            nt->_v_node[in]->eqn_index_ = i; // 1-indexed
             if (nt->_v_node[in]->extnode) {
                 i += nlayer;
             }
         }
         for (in = 0; in < nt->end; ++in) {
             int ie, k;
-            Node *nd, *pnd;
+            Node *nd;
             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) {
+                nde->eqn_index_ = i; // 0-indexed
                 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);
                 }
             }
-            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);
-                    }
-            } else { /* not needed if index starts at 1 */
-                nd->_a_matelm = nullptr;
-                nd->_b_matelm = nullptr;
-            }
         }
         nrndae_alloc();
     } else {