diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp
index 040646567a0..717994879dc 100644
--- a/src/nrniv/nonlinz.cpp
+++ b/src/nrniv/nonlinz.cpp
@@ -6,16 +6,15 @@
 #include "nrniv_mf.h"
 #include "nrnoc2iv.h"
 #include "nrnmpi.h"
-#include "cspmatrix.h"
 #include "membfunc.h"
 
 #include <Eigen/Sparse>
+#include <complex>
 
 extern void v_setup_vectors();
 extern int nrndae_extra_eqn_count();
 extern Symlist* hoc_built_in_symlist;
 extern void (*nrnthread_v_transfer_)(NrnThread*);
-extern spREAL* spGetElement(char*, int, int);
 
 extern void pargap_jacobi_rhs(double*, double*);
 extern void pargap_jacobi_setup(int mode);
@@ -36,7 +35,7 @@ class NonLinImpRep {
     int n_v_, n_ext_, n_lin_, n_ode_, neq_v_, neq_;
     std::vector<neuron::container::data_handle<double>> pv_, pvdot_;
     int* v_index_;
-    double* rv_;
+    std::vector<std::complex<double>> rv_;
     double* jv_;
     std::vector<std::complex<double>*> diag_;
     double* deltavec_;  // just like cvode.atol*cvode.atolscale for ode's
@@ -66,9 +65,7 @@ double NonLinImp::transfer_amp(int curloc, int vloc) {
     if (curloc != rep_->iloc_) {
         solve(curloc);
     }
-    double x = rep_->rv_[vloc];
-    double y = rep_->jv_[vloc];
-    return sqrt(x * x + y * y);
+    return std::abs(rep_->rv_[vloc]);
 }
 double NonLinImp::input_amp(int curloc) {
     if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_) {
@@ -80,9 +77,7 @@ double NonLinImp::input_amp(int curloc) {
     if (curloc < 0) {
         return 0.0;
     }
-    double x = rep_->rv_[curloc];
-    double y = rep_->jv_[curloc];
-    return sqrt(x * x + y * y);
+    return std::abs(rep_->rv_[curloc]);
 }
 double NonLinImp::transfer_phase(int curloc, int vloc) {
     if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_ && curloc != rep_->iloc_) {
@@ -92,9 +87,7 @@ double NonLinImp::transfer_phase(int curloc, int vloc) {
     if (curloc != rep_->iloc_) {
         solve(curloc);
     }
-    double x = rep_->rv_[vloc];
-    double y = rep_->jv_[vloc];
-    return atan2(y, x);
+    return std::arg(rep_->rv_[vloc]);
 }
 double NonLinImp::input_phase(int curloc) {
     if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_) {
@@ -106,9 +99,7 @@ double NonLinImp::input_phase(int curloc) {
     if (curloc < 0) {
         return 0.0;
     }
-    double x = rep_->rv_[curloc];
-    double y = rep_->jv_[curloc];
-    return atan2(y, x);
+    return std::arg(rep_->rv_[curloc]);
 }
 double NonLinImp::ratio_amp(int clmploc, int vloc) {
     if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_) {
@@ -120,15 +111,7 @@ double NonLinImp::ratio_amp(int clmploc, int vloc) {
     if (clmploc != rep_->iloc_) {
         solve(clmploc);
     }
-    double ax, bx, cx, ay, by, cy, bb;
-    ax = rep_->rv_[vloc];
-    ay = rep_->jv_[vloc];
-    bx = rep_->rv_[clmploc];
-    by = rep_->jv_[clmploc];
-    bb = bx * bx + by * by;
-    cx = (ax * bx + ay * by) / bb;
-    cy = (ay * bx - ax * by) / bb;
-    return sqrt(cx * cx + cy * cy);
+    return std::abs(rep_->rv_[vloc] * std::conj(rep_->rv_[clmploc]) / std::norm(rep_->rv_[clmploc]));
 }
 void NonLinImp::compute(double omega, double deltafac, int maxiter) {
     v_setup_vectors();
@@ -175,8 +158,7 @@ int NonLinImp::solve(int curloc) {
         int i;
         rep_->iloc_ = curloc;
         for (i = 0; i < rep_->neq_; ++i) {
-            rep_->rv_[i] = 0;
-            rep_->jv_[i] = 0;
+            rep_->rv_[i] = std::complex<double>(0, 0);
         }
         if (curloc >= 0) {
             rep_->rv_[curloc] = 1.e2 / NODEAREA(_nt->_v_node[curloc]);
@@ -184,20 +166,11 @@ int NonLinImp::solve(int curloc) {
         if (nrnthread_v_transfer_) {
             rval = rep_->gapsolve();
         } else {
-            std::vector<std::complex<double>> rjv{};
-            rjv.reserve(rep_->neq_);
-            for (size_t i = 0; i < rep_->neq_; ++i) {
-                rjv.emplace_back(rep_->rv_[i], rep_->jv_[i]);
-            }
             auto lu = Eigen::SparseLU<decltype(rep_->m_)>(rep_->m_);
             lu.analyzePattern(rep_->m_);
             lu.factorize(rep_->m_);
-            auto rjv_ = Eigen::Map<Eigen::Vector<std::complex<double>, Eigen::Dynamic>>(rjv.data(), rjv.size());
-            rjv_ = lu.solve(rjv_);
-            for (size_t i = 0; i < rep_->neq_; ++i) {
-                rep_->rv_[i] = rjv[i].real();
-                rep_->jv_[i] = rjv[i].imag();
-            }
+            auto rv = Eigen::Map<Eigen::Vector<std::complex<double>, Eigen::Dynamic>>(rep_->rv_.data(), rep_->rv_.size());
+            rv = lu.solve(rv);
         }
     }
     return rval;
@@ -248,10 +221,7 @@ NonLinImpRep::NonLinImpRep() {
     pv_.resize(neq_);
     pvdot_.resize(neq_);
     v_index_ = new int[n_v_];
-    rv_ = new double[neq_ + 1];
-    rv_ += 1;
-    jv_ = new double[neq_ + 1];
-    jv_ += 1;
+    rv_.resize(neq_);
     diag_.resize(neq_);
     deltavec_ = new double[neq_];
 
@@ -273,8 +243,6 @@ NonLinImpRep::NonLinImpRep() {
 
 NonLinImpRep::~NonLinImpRep() {
     delete[] v_index_;
-    delete[](rv_ - 1);
-    delete[](jv_ - 1);
     delete[] deltavec_;
 }
 
@@ -381,8 +349,6 @@ void NonLinImpRep::dids() {
             nrn_ode_count_t s = memb_func[i].ode_count;
             int cnt = (*s)(i);
             if (memb_func[i].current) {
-                double* x1 = rv_;  // use as temporary storage
-                double* x2 = jv_;
                 for (in = 0; in < ml->nodecount; ++in) {
                     Node* nd = ml->nodelist[in];
                     // zero rhs
@@ -390,18 +356,18 @@ void NonLinImpRep::dids() {
                     // compute rhs
                     current(i, ml, in);
                     // save rhs
-                    x2[in] = NODERHS(nd);
+                    rv_[in].imag() = NODERHS(nd);
                     // each state incremented separately and restored
                     for (iis = 0; iis < cnt; ++iis) {
                         is = ieq + in * cnt + iis;
                         // save s
-                        x1[is] = *pv_[is];
+                        rv_[is].real() = *pv_[is];
                         // increment s and zero rhs
                         *pv_[is] += deltavec_[is];
                         NODERHS(nd) = 0;
                         current(i, ml, in);
-                        *pv_[is] = x1[is];  // restore s
-                        double g = (NODERHS(nd) - x2[in]) / deltavec_[is];
+                        *pv_[is] = rv_[is].real();  // restore s
+                        double g = (NODERHS(nd) - rv_[in].imag()) / deltavec_[is];
                         if (g != 0.) {
                             m_.coeffRef(v_index_[nd->v_node_index], is + 1) = -g;
                         }
@@ -428,22 +394,20 @@ void NonLinImpRep::dsdv() {
             nrn_ode_count_t s = memb_func[i].ode_count;
             int cnt = (*s)(i);
             if (memb_func[i].current) {
-                double* x1 = rv_;  // use as temporary storage
-                double* x2 = jv_;
                 // zero rhs, save v
                 for (in = 0; in < ml->nodecount; ++in) {
                     Node* nd = ml->nodelist[in];
                     for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) {
                         *pvdot_[is] = 0.;
                     }
-                    x1[in] = NODEV(nd);
+                    rv_[in].real() = NODEV(nd);
                 }
                 // increment v only once in case there are multiple
                 // point processes at the same location
                 for (in = 0; in < ml->nodecount; ++in) {
                     Node* nd = ml->nodelist[in];
                     auto const v = nd->v();
-                    if (x1[in] == v) {
+                    if (rv_[in].real() == v) {
                         nd->v() = v + delta_;
                     }
                 }
@@ -453,10 +417,10 @@ void NonLinImpRep::dsdv() {
                 for (in = 0; in < ml->nodecount; ++in) {
                     Node* nd = ml->nodelist[in];
                     for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) {
-                        x2[is] = *pvdot_[is];
+                        rv_[is].imag() = *pvdot_[is];
                         *pvdot_[is] = 0;
                     }
-                    nd->v() = x1[in];
+                    nd->v() = rv_[in].real();
                 }
                 // compute the rhs(v)
                 ode(i, ml);
@@ -464,7 +428,7 @@ void NonLinImpRep::dsdv() {
                 for (in = 0; in < ml->nodecount; ++in) {
                     Node* nd = ml->nodelist[in];
                     for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) {
-                        double ds = (x2[is] - *pvdot_[is]) / delta_;
+                        double ds = (rv_[is].imag() - *pvdot_[is]) / delta_;
                         if (ds != 0.) {
                             m_.coeffRef(is + 1, v_index_[nd->v_node_index]) = -ds;
                         }
@@ -491,13 +455,11 @@ void NonLinImpRep::dsds() {
             int nc = ml->nodecount;
             nrn_ode_count_t s = memb_func[i].ode_count;
             int cnt = (*s)(i);
-            double* x1 = rv_;  // use as temporary storage
-            double* x2 = jv_;
             // zero rhs, save s
             for (in = 0; in < ml->nodecount; ++in) {
                 for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) {
                     *pvdot_[is] = 0.;
-                    x1[is] = *pv_[is];
+                    rv_[is].real() = *pv_[is];
                 }
             }
             // compute rhs. this is the rhs(s)
@@ -505,7 +467,7 @@ void NonLinImpRep::dsds() {
             // save rhs
             for (in = 0; in < ml->nodecount; ++in) {
                 for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) {
-                    x2[is] = *pvdot_[is];
+                    rv_[is].imag() = *pvdot_[is];
                 }
             }
             // iterate over the states
@@ -525,11 +487,11 @@ void NonLinImpRep::dsds() {
                     Node* nd = ml->nodelist[in];
                     ks = ieq + in * cnt + kks;
                     for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) {
-                        double ds = (*pvdot_[is] - x2[is]) / deltavec_[is];
+                        double ds = (*pvdot_[is] - rv_[is].imag()) / deltavec_[is];
                         if (ds != 0.) {
                             m_.coeffRef(is, ks) = -ds;
                         }
-                        *pv_[ks] = x1[ks];
+                        *pv_[ks] = rv_[ks].real();
                     }
                 }
                 // perhaps not necessary but ensures the last computation is with
@@ -580,22 +542,9 @@ int NonLinImpRep::gapsolve() {
 
     pargap_jacobi_setup(0);
 
-    double *rx, *jx, *rx1, *jx1, *rb, *jb;
-    if (neq_) {
-        rx = new double[neq_];
-        jx = new double[neq_];
-        rx1 = new double[neq_];
-        jx1 = new double[neq_];
-        rb = new double[neq_];
-        jb = new double[neq_];
-    }
-
-    // initialize for first iteration
-    for (int i = 0; i < neq_; ++i) {
-        rx[i] = jx[i] = 0.0;
-        rb[i] = rv_[i];
-        jb[i] = jv_[i];
-    }
+    std::vector<std::complex<double>> rx(neq_);
+    std::vector<std::complex<double>> rx1(neq_);
+    std::vector<std::complex<double>> rb(rv_.begin(), rv_.end());
 
     // iterate till change in x is small
     double tol = 1e-9;
@@ -614,19 +563,17 @@ int NonLinImpRep::gapsolve() {
             auto lu = Eigen::SparseLU<decltype(m_)>(m_);
             lu.analyzePattern(m_);
             lu.factorize(m_);
-            auto rjv_ = Eigen::Map<Eigen::Vector<std::complex<double>, Eigen::Dynamic>>(rjv.data(), rjv.size());
-            rjv_ = lu.solve(rjv_);
-            for (size_t i = 0; i < neq_; ++i) {
-                jb[i] = rjv[i].real();
-                jx1[i] = rjv[i].imag();
-            }
+            auto rb_ = Eigen::Map<Eigen::Vector<std::complex<double>, Eigen::Dynamic>>(rb.data(), rb.size());
+            auto rx1_ = Eigen::Map<Eigen::Vector<std::complex<double>, Eigen::Dynamic>>(rx1.data(), rx1.size());
+            rx1_ = lu.solve(rb_);
         }
 
         // if any change in x > tol, then do another iteration.
         success = 1;
         delta = 0.0;
         for (int i = 0; i < neq_; ++i) {
-            double err = fabs(rx1[i] - rx[i]) + fabs(jx1[i] - jx[i]);
+            auto diff = rx1[i] - rx[i];
+            double err = std::abs(diff.real()) + std::abs(diff.imag());
             if (err > tol) {
                 success = 0;
             }
@@ -640,20 +587,13 @@ int NonLinImpRep::gapsolve() {
         }
 #endif
         if (success) {
-            for (int i = 0; i < neq_; ++i) {
-                rv_[i] = rx1[i];
-                jv_[i] = jx1[i];
-            }
+            rv_ = rx1;
             break;
         }
 
         // setup for next iteration
-        for (int i = 0; i < neq_; ++i) {
-            rx[i] = rx1[i];
-            jx[i] = jx1[i];
-            rb[i] = rv_[i];
-            jb[i] = jv_[i];
-        }
+        rx = rx1;
+        rb = rv_;
         pargap_jacobi_rhs(rb, rx);
         pargap_jacobi_rhs(jb, jx);
     }