diff --git a/src/ivoc/matrix.cpp b/src/ivoc/matrix.cpp
index 8fa098d61d..ddc8d41cfa 100644
--- a/src/ivoc/matrix.cpp
+++ b/src/ivoc/matrix.cpp
@@ -64,13 +64,10 @@ static double m_ncol(void* v) {
 
 static double m_setval(void* v) {
     Matrix* m = (Matrix*) v;
-    int i, j;
-    double val, *pval;
-    i = (int) chkarg(1, 0, m->nrow() - 1);
-    j = (int) chkarg(2, 0, m->ncol() - 1);
-    val = *getarg(3);
-    pval = m->mep(i, j);
-    *pval = val;
+    int i = (int) chkarg(1, 0, m->nrow() - 1);
+    int j = (int) chkarg(2, 0, m->ncol() - 1);
+    double val = *getarg(3);
+    m->coeff(i, j) = val;
     return val;
 }
 
@@ -171,7 +168,7 @@ static double m_scanf(void* v) {
     m->resize(nrow, ncol);
     for (i = 0; i < nrow; ++i)
         for (j = 0; j < ncol; ++j) {
-            *(m->mep(i, j)) = hoc_scan(f);
+            m->coeff(i, j) = hoc_scan(f);
         }
     return 0.;
 }
@@ -602,7 +599,7 @@ static Object** m_set(void* v) {
     int k;
     for (k = 0, i = 0; i < nrow; ++i) {
         for (j = 0; j < ncol; ++j) {
-            *(m->mep(i, j)) = *getarg(++k);
+            m->coeff(i, j) = *getarg(++k);
         }
     }
     return temp_objvar(m);
@@ -641,7 +638,7 @@ static Object** m_from_vector(void* v) {
     double* ve = vector_vec(vout);
     for (j = 0; j < ncol; ++j)
         for (i = 0; i < nrow; ++i) {
-            *(m->mep(i, j)) = ve[k++];
+            m->coeff(i, j) = ve[k++];
         }
     return temp_objvar(m);
 }
diff --git a/src/ivoc/ocmatrix.cpp b/src/ivoc/ocmatrix.cpp
index 7afa87498a..593bdef793 100644
--- a/src/ivoc/ocmatrix.cpp
+++ b/src/ivoc/ocmatrix.cpp
@@ -34,21 +34,20 @@ OcMatrix* OcMatrix::instance(int nrow, int ncol, int type) {
     }
 }
 
-void OcMatrix::unimp() {
-    hoc_execerror("Matrix method not implemented for this type matrix", 0);
+void OcMatrix::unimp() const {
+    hoc_execerror("Matrix method not implemented for this type matrix", nullptr);
 }
 
-void OcMatrix::nonzeros(std::vector<int>& m, std::vector<int>& n) {
-    m.clear();
-    n.clear();
+std::vector<std::pair<int, int>> OcMatrix::nonzeros() const {
+    std::vector<std::pair<int, int>> nzs;
     for (int i = 0; i < nrow(); i++) {
         for (int j = 0; j < ncol(); j++) {
-            if (getval(i, j) != 0) {
-                m.push_back(i);
-                n.push_back(j);
+            if (getval(i, j) != 0.) {
+                nzs.emplace_back(i, j);
             }
         }
     }
+    return nzs;
 }
 
 OcFullMatrix* OcMatrix::full() {
@@ -65,16 +64,16 @@ OcFullMatrix::OcFullMatrix(int nrow, int ncol)
     m_.setZero();
 }
 
-double* OcFullMatrix::mep(int i, int j) {
-    return &m_(i, j);
+double& OcFullMatrix::coeff(int i, int j) {
+    return m_(i, j);
 }
-double OcFullMatrix::getval(int i, int j) {
+double OcFullMatrix::getval(int i, int j) const {
     return m_(i, j);
 }
-int OcFullMatrix::nrow() {
+int OcFullMatrix::nrow() const {
     return m_.rows();
 }
-int OcFullMatrix::ncol() {
+int OcFullMatrix::ncol() const {
     return m_.cols();
 }
 
@@ -84,29 +83,29 @@ void OcFullMatrix::resize(int i, int j) {
     m_.conservativeResizeLike(v);
 }
 
-void OcFullMatrix::mulv(Vect* vin, Vect* vout) {
+void OcFullMatrix::mulv(Vect* vin, Vect* vout) const {
     auto v1 = Vect2VEC(vin);
     auto v2 = Vect2VEC(vout);
     v2 = m_ * v1;
 }
 
-void OcFullMatrix::mulm(Matrix* in, Matrix* out) {
+void OcFullMatrix::mulm(Matrix* in, Matrix* out) const {
     out->full()->m_ = m_ * in->full()->m_;
 }
 
-void OcFullMatrix::muls(double s, Matrix* out) {
+void OcFullMatrix::muls(double s, Matrix* out) const {
     out->full()->m_ = s * m_;
 }
 
-void OcFullMatrix::add(Matrix* in, Matrix* out) {
+void OcFullMatrix::add(Matrix* in, Matrix* out) const {
     out->full()->m_ = m_ + in->full()->m_;
 }
 
-void OcFullMatrix::copy(Matrix* out) {
+void OcFullMatrix::copy(Matrix* out) const {
     out->full()->m_ = m_;
 }
 
-void OcFullMatrix::bcopy(Matrix* out, int i0, int j0, int n0, int m0, int i1, int j1) {
+void OcFullMatrix::bcopy(Matrix* out, int i0, int j0, int n0, int m0, int i1, int j1) const {
     out->full()->m_.block(i1, j1, n0, m0) = m_.block(i0, j0, n0, m0);
 }
 
@@ -119,14 +118,14 @@ void OcFullMatrix::transpose(Matrix* out) {
 }
 
 // As only symmetric matrix are accepted, eigenvalues are not complex
-void OcFullMatrix::symmeigen(Matrix* mout, Vect* vout) {
+void OcFullMatrix::symmeigen(Matrix* mout, Vect* vout) const {
     auto v1 = Vect2VEC(vout);
     Eigen::EigenSolver<Eigen::MatrixXd> es(m_);
     v1 = es.eigenvalues().real();
     mout->full()->m_ = es.eigenvectors().real();
 }
 
-void OcFullMatrix::svd1(Matrix* u, Matrix* v, Vect* d) {
+void OcFullMatrix::svd1(Matrix* u, Matrix* v, Vect* d) const {
     auto v1 = Vect2VEC(d);
     Eigen::JacobiSVD<Eigen::MatrixXd> svd(m_, Eigen::ComputeFullU | Eigen::ComputeFullV);
     v1 = svd.singularValues();
@@ -138,17 +137,17 @@ void OcFullMatrix::svd1(Matrix* u, Matrix* v, Vect* d) {
     }
 }
 
-void OcFullMatrix::getrow(int k, Vect* out) {
+void OcFullMatrix::getrow(int k, Vect* out) const {
     auto v1 = Vect2VEC(out);
     v1 = m_.row(k);
 }
 
-void OcFullMatrix::getcol(int k, Vect* out) {
+void OcFullMatrix::getcol(int k, Vect* out) const {
     auto v1 = Vect2VEC(out);
     v1 = m_.col(k);
 }
 
-void OcFullMatrix::getdiag(int k, Vect* out) {
+void OcFullMatrix::getdiag(int k, Vect* out) const {
     auto vout = m_.diagonal(k);
     if (k >= 0) {
         for (int i = 0, j = k; i < nrow() && j < ncol(); ++i, ++j) {
@@ -205,15 +204,15 @@ void OcFullMatrix::ident() {
     m_.setIdentity();
 }
 
-void OcFullMatrix::exp(Matrix* out) {
+void OcFullMatrix::exp(Matrix* out) const {
     out->full()->m_ = m_.exp();
 }
 
-void OcFullMatrix::pow(int i, Matrix* out) {
+void OcFullMatrix::pow(int i, Matrix* out) const {
     out->full()->m_ = m_.pow(i).eval();
 }
 
-void OcFullMatrix::inverse(Matrix* out) {
+void OcFullMatrix::inverse(Matrix* out) const {
     out->full()->m_ = m_.inverse();
 }
 
@@ -226,7 +225,7 @@ void OcFullMatrix::solv(Vect* in, Vect* out, bool use_lu) {
     v2 = lu_->solve(v1);
 }
 
-double OcFullMatrix::det(int* e) {
+double OcFullMatrix::det(int* e) const {
     *e = 0;
     double m = m_.determinant();
     if (m) {
@@ -248,8 +247,8 @@ OcSparseMatrix::OcSparseMatrix(int nrow, int ncol)
     : OcMatrix(MSPARSE)
     , m_(nrow, ncol) {}
 
-double* OcSparseMatrix::mep(int i, int j) {
-    return &m_.coeffRef(i, j);
+double& OcSparseMatrix::coeff(int i, int j) {
+    return m_.coeffRef(i, j);
 }
 
 void OcSparseMatrix::zero() {
@@ -260,19 +259,19 @@ void OcSparseMatrix::zero() {
     }
 }
 
-double OcSparseMatrix::getval(int i, int j) {
+double OcSparseMatrix::getval(int i, int j) const {
     return m_.coeff(i, j);
 }
 
-int OcSparseMatrix::nrow() {
+int OcSparseMatrix::nrow() const {
     return m_.rows();
 }
 
-int OcSparseMatrix::ncol() {
+int OcSparseMatrix::ncol() const {
     return m_.cols();
 }
 
-void OcSparseMatrix::mulv(Vect* vin, Vect* vout) {
+void OcSparseMatrix::mulv(Vect* vin, Vect* vout) const {
     auto v1 = Vect2VEC(vin);
     auto v2 = Vect2VEC(vout);
     v2 = m_ * v1;
@@ -330,7 +329,7 @@ void OcSparseMatrix::setcol(int k, double in) {
     }
 }
 
-void OcSparseMatrix::ident(void) {
+void OcSparseMatrix::ident() {
     m_.setIdentity();
 }
 
@@ -348,7 +347,7 @@ void OcSparseMatrix::setdiag(int k, double in) {
     }
 }
 
-int OcSparseMatrix::sprowlen(int i) {
+int OcSparseMatrix::sprowlen(int i) const {
     int acc = 0;
     for (decltype(m_)::InnerIterator it(m_, i); it; ++it) {
         acc += 1;
@@ -356,7 +355,7 @@ int OcSparseMatrix::sprowlen(int i) {
     return acc;
 }
 
-double OcSparseMatrix::spgetrowval(int i, int jindx, int* j) {
+double OcSparseMatrix::spgetrowval(int i, int jindx, int* j) const {
     int acc = 0;
     for (decltype(m_)::InnerIterator it(m_, i); it; ++it) {
         if (acc == jindx) {
@@ -368,13 +367,13 @@ double OcSparseMatrix::spgetrowval(int i, int jindx, int* j) {
     return 0;
 }
 
-void OcSparseMatrix::nonzeros(std::vector<int>& m, std::vector<int>& n) {
-    m.clear();
-    n.clear();
+std::vector<std::pair<int, int>> OcSparseMatrix::nonzeros() const {
+    std::vector<std::pair<int, int>> nzs;
+    nzs.reserve(m_.nonZeros());
     for (int k = 0; k < m_.outerSize(); ++k) {
         for (decltype(m_)::InnerIterator it(m_, k); it; ++it) {
-            m.push_back(it.row());
-            n.push_back(it.col());
+            nzs.emplace_back(it.row(), it.col());
         }
     }
+    return nzs;
 }
diff --git a/src/ivoc/ocmatrix.h b/src/ivoc/ocmatrix.h
index 15960fe51f..9e874c86a3 100644
--- a/src/ivoc/ocmatrix.h
+++ b/src/ivoc/ocmatrix.h
@@ -1,6 +1,7 @@
 #pragma once
 
 #include <memory>
+#include <utility>
 #include <vector>
 
 #include <Eigen/Eigen>
@@ -20,23 +21,35 @@ class OcMatrix {
     static OcMatrix* instance(int nrow, int ncol, int type = MFULL);
     virtual ~OcMatrix() = default;
 
-    virtual double* mep(int i, int j) {
+    // This function is deprecated and should not be used!
+    // mep stands for 'matrix element pointer'
+    inline double* mep(int i, int j) {
+        return &coeff(i, j);
+    }
+
+    inline double operator()(int i, int j) const {
+        return getval(i, j);
+    };
+
+    virtual double& coeff(int i, int j) {
+        static double zero = 0.0;
         unimp();
-        return nullptr;
-    }  // matrix element pointer
+        return zero;
+    }
+
     inline double& operator()(int i, int j) {
-        return *mep(i, j);
+        return coeff(i, j);
     };
 
-    virtual double getval(int i, int j) {
+    virtual double getval(int i, int j) const {
         unimp();
         return 0.;
     }
-    virtual int nrow() {
+    virtual int nrow() const {
         unimp();
         return 0;
     }
-    virtual int ncol() {
+    virtual int ncol() const {
         unimp();
         return 0;
     }
@@ -44,32 +57,32 @@ class OcMatrix {
         unimp();
     }
 
-    virtual void nonzeros(std::vector<int>& m, std::vector<int>& n);
+    virtual std::vector<std::pair<int, int>> nonzeros() const;
 
     OcFullMatrix* full();
 
-    inline void mulv(Vect& in, Vect& out) {
+    inline void mulv(Vect& in, Vect& out) const {
         mulv(&in, &out);
     };
-    virtual void mulv(Vect* in, Vect* out) {
+    virtual void mulv(Vect* in, Vect* out) const {
         unimp();
     }
-    virtual void mulm(Matrix* in, Matrix* out) {
+    virtual void mulm(Matrix* in, Matrix* out) const {
         unimp();
     }
-    virtual void muls(double, Matrix* out) {
+    virtual void muls(double, Matrix* out) const {
         unimp();
     }
-    virtual void add(Matrix*, Matrix* out) {
+    virtual void add(Matrix*, Matrix* out) const {
         unimp();
     }
-    virtual void getrow(int, Vect* out) {
+    virtual void getrow(int, Vect* out) const {
         unimp();
     }
-    virtual void getcol(int, Vect* out) {
+    virtual void getcol(int, Vect* out) const {
         unimp();
     }
-    virtual void getdiag(int, Vect* out) {
+    virtual void getdiag(int, Vect* out) const {
         unimp();
     }
     virtual void setrow(int, Vect* in) {
@@ -96,47 +109,47 @@ class OcMatrix {
     virtual void ident() {
         unimp();
     }
-    virtual void exp(Matrix* out) {
+    virtual void exp(Matrix* out) const {
         unimp();
     }
-    virtual void pow(int, Matrix* out) {
+    virtual void pow(int, Matrix* out) const {
         unimp();
     }
-    virtual void inverse(Matrix* out) {
+    virtual void inverse(Matrix* out) const {
         unimp();
     }
     virtual void solv(Vect* vin, Vect* vout, bool use_lu) {
         unimp();
     }
-    virtual void copy(Matrix* out) {
+    virtual void copy(Matrix* out) const {
         unimp();
     }
-    virtual void bcopy(Matrix* mout, int i0, int j0, int n0, int m0, int i1, int j1) {
+    virtual void bcopy(Matrix* mout, int i0, int j0, int n0, int m0, int i1, int j1) const {
         unimp();
     }
     virtual void transpose(Matrix* out) {
         unimp();
     }
-    virtual void symmeigen(Matrix* mout, Vect* vout) {
+    virtual void symmeigen(Matrix* mout, Vect* vout) const {
         unimp();
     }
-    virtual void svd1(Matrix* u, Matrix* v, Vect* d) {
+    virtual void svd1(Matrix* u, Matrix* v, Vect* d) const {
         unimp();
     }
-    virtual double det(int* e) {
+    virtual double det(int* e) const {
         unimp();
         return 0.0;
     }
-    virtual int sprowlen(int) {
+    virtual int sprowlen(int) const {
         unimp();
         return 0;
     }
-    virtual double spgetrowval(int i, int jindx, int* j) {
+    virtual double spgetrowval(int i, int jindx, int* j) const {
         unimp();
         return 0.;
     }
 
-    void unimp();
+    void unimp() const;
 
   protected:
     OcMatrix(int type);
@@ -155,19 +168,19 @@ class OcFullMatrix final: public OcMatrix {  // type 1
     OcFullMatrix(int, int);
     ~OcFullMatrix() override = default;
 
-    double* mep(int, int) override;
-    double getval(int i, int j) override;
-    int nrow() override;
-    int ncol() override;
+    double& coeff(int, int) override;
+    double getval(int i, int j) const override;
+    int nrow() const override;
+    int ncol() const override;
     void resize(int, int) override;
 
-    void mulv(Vect* in, Vect* out) override;
-    void mulm(Matrix* in, Matrix* out) override;
-    void muls(double, Matrix* out) override;
-    void add(Matrix*, Matrix* out) override;
-    void getrow(int, Vect* out) override;
-    void getcol(int, Vect* out) override;
-    void getdiag(int, Vect* out) override;
+    void mulv(Vect* in, Vect* out) const override;
+    void mulm(Matrix* in, Matrix* out) const override;
+    void muls(double, Matrix* out) const override;
+    void add(Matrix*, Matrix* out) const override;
+    void getrow(int, Vect* out) const override;
+    void getcol(int, Vect* out) const override;
+    void getdiag(int, Vect* out) const override;
     void setrow(int, Vect* in) override;
     void setcol(int, Vect* in) override;
     void setdiag(int, Vect* in) override;
@@ -176,16 +189,16 @@ class OcFullMatrix final: public OcMatrix {  // type 1
     void setdiag(int, double in) override;
     void zero() override;
     void ident() override;
-    void exp(Matrix* out) override;
-    void pow(int, Matrix* out) override;
-    void inverse(Matrix* out) override;
+    void exp(Matrix* out) const override;
+    void pow(int, Matrix* out) const override;
+    void inverse(Matrix* out) const override;
     void solv(Vect* vin, Vect* vout, bool use_lu) override;
-    void copy(Matrix* out) override;
-    void bcopy(Matrix* mout, int i0, int j0, int n0, int m0, int i1, int j1) override;
+    void copy(Matrix* out) const override;
+    void bcopy(Matrix* mout, int i0, int j0, int n0, int m0, int i1, int j1) const override;
     void transpose(Matrix* out) override;
-    void symmeigen(Matrix* mout, Vect* vout) override;
-    void svd1(Matrix* u, Matrix* v, Vect* d) override;
-    double det(int* exponent) override;
+    void symmeigen(Matrix* mout, Vect* vout) const override;
+    void svd1(Matrix* u, Matrix* v, Vect* d) const override;
+    double det(int* exponent) const override;
 
   private:
     Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> m_{};
@@ -197,12 +210,12 @@ class OcSparseMatrix final: public OcMatrix {  // type 2
     OcSparseMatrix(int, int);
     ~OcSparseMatrix() override = default;
 
-    double* mep(int, int) override;
-    int nrow() override;
-    int ncol() override;
-    double getval(int, int) override;
-    void ident(void) override;
-    void mulv(Vect* in, Vect* out) override;
+    double& coeff(int, int) override;
+    int nrow() const override;
+    int ncol() const override;
+    double getval(int, int) const override;
+    void ident() override;
+    void mulv(Vect* in, Vect* out) const override;
     void solv(Vect* vin, Vect* vout, bool use_lu) override;
 
     void setrow(int, Vect* in) override;
@@ -212,10 +225,10 @@ class OcSparseMatrix final: public OcMatrix {  // type 2
     void setcol(int, double in) override;
     void setdiag(int, double in) override;
 
-    void nonzeros(std::vector<int>& m, std::vector<int>& n) override;
+    std::vector<std::pair<int, int>> nonzeros() const override;
 
-    int sprowlen(int) override;  // how many elements in row
-    double spgetrowval(int i, int jindx, int* j) override;
+    int sprowlen(int) const override;  // how many elements in row
+    double spgetrowval(int i, int jindx, int* j) const override;
 
     void zero() override;
 
diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp
index f07cfbed2b..658e66c5c0 100644
--- a/src/nrniv/matrixmap.cpp
+++ b/src/nrniv/matrixmap.cpp
@@ -19,7 +19,8 @@ void MatrixMap::mmfree() {
 
 void MatrixMap::add(double fac) {
     for (int i = 0; i < plen_; ++i) {
-        *ptree_[i] += fac * (*pm_[i]);
+        auto [it, jt] = pm_[i];
+        *ptree_[i] += fac * m_(it, jt);
     }
 }
 
@@ -28,13 +29,10 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) {
     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());
-    for (int k = 0; k < nonzero_i.size(); k++) {
-        const int i = nonzero_i[k];
-        const int j = nonzero_j[k];
+    std::vector<std::pair<int, int>> nzs = m_.nonzeros();
+    pm_.resize(nzs.size());
+    ptree_.resize(nzs.size());
+    for (const auto [i, j]: nzs) {
         int it;
         if (i < nnode) {
             it = nodes[i]->eqn_index_ + layer[i];
@@ -45,7 +43,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_[plen_] = std::make_pair(i, j);
         if (j < nnode) {
             jt = nodes[j]->eqn_index_ + layer[j];
             if (layer[j] > 0 && !nodes[j]->extnode) {
diff --git a/src/nrniv/matrixmap.h b/src/nrniv/matrixmap.h
index f216dcc9c3..c20b774327 100644
--- a/src/nrniv/matrixmap.h
+++ b/src/nrniv/matrixmap.h
@@ -116,6 +116,6 @@ class MatrixMap {
 
     // the map
     int plen_ = 0;
-    std::vector<double*> pm_{};
+    std::vector<std::pair<int, int>> pm_{};
     std::vector<double*> ptree_{};
 };
diff --git a/src/nrnpython/rxd.cpp b/src/nrnpython/rxd.cpp
index 7cf61f7513..ff548a7f4a 100644
--- a/src/nrnpython/rxd.cpp
+++ b/src/nrnpython/rxd.cpp
@@ -1516,7 +1516,7 @@ void solve_reaction(ICSReactions* react,
     double pd;
     double dt = *dt_ptr;
     double dx = FLT_EPSILON;
-    auto jacobian = std::make_unique<OcFullMatrix>(N, N);
+    OcFullMatrix jacobian(N, N);
     auto b = std::make_unique<IvocVect>(N);
     auto x = std::make_unique<IvocVect>(N);
 
@@ -1655,7 +1655,7 @@ void solve_reaction(ICSReactions* react,
                             if (react->state_idx[segment][jac_i][jac_j] != SPECIES_ABSENT) {
                                 pd = (result_array_dx[jac_i][jac_j] - result_array[jac_i][jac_j]) /
                                      dx;
-                                *jacobian->mep(jac_idx, idx) = (idx == jac_idx) - dt * pd;
+                                jacobian(jac_idx, idx) = (idx == jac_idx) - dt * pd;
                                 jac_idx += 1;
                             }
                             result_array_dx[jac_i][jac_j] = 0;
@@ -1665,7 +1665,7 @@ void solve_reaction(ICSReactions* react,
                         // pd is our Jacobian approximated
                         if (react->ecs_state[segment][jac_i] != NULL) {
                             pd = (ecs_result_dx[jac_i] - ecs_result[jac_i]) / dx;
-                            *jacobian->mep(jac_idx, idx) = -dt * pd;
+                            jacobian(jac_idx, idx) = -dt * pd;
                             jac_idx += 1;
                         }
                         ecs_result_dx[jac_i] = 0;
@@ -1707,7 +1707,7 @@ void solve_reaction(ICSReactions* react,
                         // pd is our Jacobian approximated
                         if (react->state_idx[segment][jac_i][jac_j] != SPECIES_ABSENT) {
                             pd = (result_array_dx[jac_i][jac_j] - result_array[jac_i][jac_j]) / dx;
-                            *jacobian->mep(jac_idx, idx) = -dt * pd;
+                            jacobian(jac_idx, idx) = -dt * pd;
                             jac_idx += 1;
                         }
                     }
@@ -1716,10 +1716,10 @@ void solve_reaction(ICSReactions* react,
                     // pd is our Jacobian approximated
                     if (react->ecs_state[segment][jac_i] != NULL) {
                         pd = (ecs_result_dx[jac_i] - ecs_result[jac_i]) / dx;
-                        *jacobian->mep(jac_idx, idx) = (idx == jac_idx) - dt * pd;
+                        jacobian(jac_idx, idx) = (idx == jac_idx) - dt * pd;
                         jac_idx += 1;
                     } else {
-                        *jacobian->mep(idx, idx) = 1.0;
+                        jacobian(idx, idx) = 1.0;
                     }
                     // reset dx array
                     ecs_states_for_reaction_dx[i] -= dx;
@@ -1728,7 +1728,7 @@ void solve_reaction(ICSReactions* react,
             }
         }
         // solve for x, destructively
-        jacobian->solv(b.get(), x.get(), false);
+        jacobian.solv(b.get(), x.get(), false);
 
         if (bval != NULL)  // variable-step
         {
diff --git a/src/nrnpython/rxd_extracellular.cpp b/src/nrnpython/rxd_extracellular.cpp
index aeb377eed4..a26d61eeb2 100644
--- a/src/nrnpython/rxd_extracellular.cpp
+++ b/src/nrnpython/rxd_extracellular.cpp
@@ -290,7 +290,6 @@ void* ecs_do_reactions(void* dataptr) {
     double* mc_mults_array = NULL;
     double dx = FLT_EPSILON;
     double pd;
-    std::unique_ptr<OcFullMatrix> jacobian;
     std::vector<double> x{};
     std::vector<double> b{};
 
@@ -316,8 +315,7 @@ void* ecs_do_reactions(void* dataptr) {
                 if (react->num_species_involved == 0)
                     continue;
                 /*allocate data structures*/
-                jacobian = std::make_unique<OcFullMatrix>(react->num_species_involved,
-                                                          react->num_species_involved);
+                OcFullMatrix jacobian(react->num_species_involved, react->num_species_involved);
                 b.resize(react->num_species_involved);
                 x.resize(react->num_species_involved);
                 states_cache = (double*) malloc(sizeof(double) * react->num_species_involved);
@@ -358,23 +356,23 @@ void* ecs_do_reactions(void* dataptr) {
 
                             for (k = 0; k < react->num_species_involved; k++) {
                                 pd = (results_array_dx[k] - results_array[k]) / dx;
-                                *jacobian->mep(k, j) = (j == k) - dt * pd;
+                                jacobian(k, j) = (j == k) - dt * pd;
                             }
                             states_cache_dx[j] -= dx;
                         }
                         // solve for x
                         if (react->num_species_involved == 1) {
-                            react->species_states[0][i] += b[0] / jacobian->getval(0, 0);
+                            react->species_states[0][i] += b[0] / jacobian(0, 0);
                         } else {
                             // find entry in leftmost column with largest absolute value
                             // Pivot
                             for (j = 0; j < react->num_species_involved; j++) {
                                 for (k = j + 1; k < react->num_species_involved; k++) {
-                                    if (abs(jacobian->getval(j, j)) < abs(jacobian->getval(k, j))) {
+                                    if (abs(jacobian(j, j)) < abs(jacobian(k, j))) {
                                         for (n = 0; n < react->num_species_involved; n++) {
-                                            temp = jacobian->getval(j, n);
-                                            *jacobian->mep(j, n) = jacobian->getval(k, n);
-                                            *jacobian->mep(k, n) = temp;
+                                            temp = jacobian(j, n);
+                                            jacobian(j, n) = jacobian(k, n);
+                                            jacobian(k, n) = temp;
                                         }
                                     }
                                 }
@@ -382,11 +380,11 @@ void* ecs_do_reactions(void* dataptr) {
 
                             for (j = 0; j < react->num_species_involved - 1; j++) {
                                 for (k = j + 1; k < react->num_species_involved; k++) {
-                                    ge_value = jacobian->getval(k, j) / jacobian->getval(j, j);
+                                    ge_value = jacobian(k, j) / jacobian(j, j);
                                     for (n = 0; n < react->num_species_involved; n++) {
-                                        val_to_set = jacobian->getval(k, n) -
-                                                     ge_value * jacobian->getval(j, n);
-                                        *jacobian->mep(k, n) = val_to_set;
+                                        val_to_set = jacobian(k, n) -
+                                                     ge_value * jacobian(j, n);
+                                        jacobian(k, n) = val_to_set;
                                     }
                                     b[k] = b[k] - ge_value * b[j];
                                 }
@@ -396,10 +394,10 @@ void* ecs_do_reactions(void* dataptr) {
                                 x[j] = b[j];
                                 for (k = j + 1; k < react->num_species_involved; k++) {
                                     if (k != j) {
-                                        x[j] = x[j] - jacobian->getval(j, k) * x[k];
+                                        x[j] = x[j] - jacobian(j, k) * x[k];
                                     }
                                 }
-                                x[j] = x[j] / jacobian->getval(j, j);
+                                x[j] = x[j] / jacobian(j, j);
                             }
                             for (j = 0; j < react->num_species_involved; j++) {
                                 // I think this should be something like
@@ -446,8 +444,7 @@ void* ecs_do_reactions(void* dataptr) {
                 if (react->num_species_involved == 0)
                     continue;
                 /*allocate data structures*/
-                jacobian = std::make_unique<OcFullMatrix>(react->num_species_involved,
-                                                          react->num_species_involved);
+                OcFullMatrix jacobian(react->num_species_involved, react->num_species_involved);
                 b.resize(react->num_species_involved);
                 x.resize(react->num_species_involved);
                 states_cache = (double*) malloc(sizeof(double) * react->num_species_involved);
@@ -480,23 +477,23 @@ void* ecs_do_reactions(void* dataptr) {
 
                             for (k = 0; k < react->num_species_involved; k++) {
                                 pd = (results_array_dx[k] - results_array[k]) / dx;
-                                *jacobian->mep(k, j) = (j == k) - dt * pd;
+                                jacobian(k, j) = (j == k) - dt * pd;
                             }
                             states_cache_dx[j] -= dx;
                         }
                         // solve for x
                         if (react->num_species_involved == 1) {
-                            react->species_states[0][i] += b[0] / jacobian->getval(0, 0);
+                            react->species_states[0][i] += b[0] / jacobian(0, 0);
                         } else {
                             // find entry in leftmost column with largest absolute value
                             // Pivot
                             for (j = 0; j < react->num_species_involved; j++) {
                                 for (k = j + 1; k < react->num_species_involved; k++) {
-                                    if (abs(jacobian->getval(j, j)) < abs(jacobian->getval(k, j))) {
+                                    if (abs(jacobian(j, j)) < abs(jacobian(k, j))) {
                                         for (n = 0; n < react->num_species_involved; n++) {
-                                            temp = jacobian->getval(j, n);
-                                            *jacobian->mep(j, n) = jacobian->getval(k, n);
-                                            *jacobian->mep(k, n) = temp;
+                                            temp = jacobian(j, n);
+                                            jacobian(j, n) = jacobian(k, n);
+                                            jacobian(k, n) = temp;
                                         }
                                     }
                                 }
@@ -504,11 +501,11 @@ void* ecs_do_reactions(void* dataptr) {
 
                             for (j = 0; j < react->num_species_involved - 1; j++) {
                                 for (k = j + 1; k < react->num_species_involved; k++) {
-                                    ge_value = jacobian->getval(k, j) / jacobian->getval(j, j);
+                                    ge_value = jacobian(k, j) / jacobian(j, j);
                                     for (n = 0; n < react->num_species_involved; n++) {
-                                        val_to_set = jacobian->getval(k, n) -
-                                                     ge_value * jacobian->getval(j, n);
-                                        *jacobian->mep(k, n) = val_to_set;
+                                        val_to_set = jacobian(k, n) -
+                                                     ge_value * jacobian(j, n);
+                                        jacobian(k, n) = val_to_set;
                                     }
                                     b[k] = b[k] - ge_value * b[j];
                                 }
@@ -518,10 +515,10 @@ void* ecs_do_reactions(void* dataptr) {
                                 x[j] = b[j];
                                 for (k = j + 1; k < react->num_species_involved; k++) {
                                     if (k != j) {
-                                        x[j] = x[j] - jacobian->getval(j, k) * x[k];
+                                        x[j] = x[j] - jacobian(j, k) * x[k];
                                     }
                                 }
-                                x[j] = x[j] / jacobian->getval(j, j);
+                                x[j] = x[j] / jacobian(j, j);
                             }
                             for (j = 0; j < react->num_species_involved; j++) {
                                 // I think this should be something like
diff --git a/test/unit_tests/matrix.cpp b/test/unit_tests/matrix.cpp
index 27386501ca..115dbfe120 100644
--- a/test/unit_tests/matrix.cpp
+++ b/test/unit_tests/matrix.cpp
@@ -94,6 +94,11 @@ SCENARIO("A Matrix", "[neuron_ivoc][OcMatrix]") {
             m.setrow(3, 2.0);
             REQUIRE(compareMatrix(m, {{3., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}, {2., 2., 2.}}));
         }
+        {
+            std::vector<std::pair<int, int>> nzs = m.nonzeros();
+            std::vector<std::pair<int, int>> res = {{0, 0}, {1, 1} , {2, 2}, {3, 0}, {3, 1}, {3, 2}};
+            REQUIRE(nzs == res);
+        }
         {
             std::vector<int> x, y;
             m.nonzeros(x, y);