diff --git a/elec-prop-test/make_figures_3D.py b/elec-prop-test/make_figures_3D.py
new file mode 100644
index 0000000..d7378b5
--- /dev/null
+++ b/elec-prop-test/make_figures_3D.py
@@ -0,0 +1,391 @@
+import matplotlib as mpl
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+import sys
+
+from fenics import * 
+import string
+
+from knpemidg import pcws_constant_project
+from knpemidg import interface_normal, plus, minus
+
+JUMP = lambda f, n: minus(f, n) - plus(f, n)
+
+# set font & text parameters
+font = {'family' : 'serif',
+        'weight' : 'bold',
+        'size'   : 13}
+
+plt.rc('font', **font)
+plt.rc('text', usetex=True)
+mpl.rcParams['image.cmap'] = 'jet'
+
+path = 'results/data/'
+
+c_1 = '#%02x%02x%02x' % (255, 0, 255)  # pink
+c_2 = '#%02x%02x%02x' % (54, 92, 141)  # blue
+c_3 = '#%02x%02x%02x' % (255, 165, 0)  # orange
+phi_M_c = '#%02x%02x%02x' % (54, 92, 141)
+
+def get_time_series(dt, T, fname, x, y, z):
+    # read data file
+    hdf5file = HDF5File(MPI.comm_world, fname, "r")
+
+    mesh = Mesh()
+    subdomains = MeshFunction("size_t", mesh, 2)
+    surfaces = MeshFunction("size_t", mesh, 1)
+    hdf5file.read(mesh, '/mesh', False)
+    mesh.coordinates()[:] *= 1e6
+    hdf5file.read(subdomains, '/subdomains')
+    hdf5file.read(surfaces, '/surfaces')
+
+    P1 = FiniteElement('CG', mesh.ufl_cell(), 1)
+    W = FunctionSpace(mesh, MixedElement(2*[P1]))
+    V = FunctionSpace(mesh, P1)
+
+    u = Function(W)
+    v = Function(V)
+    w = Function(V)
+
+    f_Na = Function(V)
+    f_K = Function(V)
+    f_Cl = Function(V)
+    f_phi = Function(V)
+
+    Na = []
+    K = []
+    Cl = []
+    phi = []
+
+    for n in range(1, int(T/dt)):
+
+            # read file
+            hdf5file.read(u, "/concentrations/vector_" + str(n))
+
+            # K concentrations
+            assign(f_K, u.sub(0))
+            K.append(f_K(x, y, z))
+
+            # Cl concentrations
+            assign(f_Cl, u.sub(1))
+            Cl.append(f_Cl(x, y, z))
+
+            # Na concentrations
+            hdf5file.read(v, "/elim_concentration/vector_" + str(n))
+            assign(f_Na, v)
+            Na.append(f_Na(x, y, z))
+
+            # potential
+            hdf5file.read(w, "/potential/vector_" + str(n))
+            assign(f_phi, w)
+            phi.append(1.0e3*f_phi(x, y, z))
+
+    return Na, K, Cl, phi
+
+def get_time_series_membrane(dt, T, fname, x_, y_, z_):
+    # read data file
+    hdf5file = HDF5File(MPI.comm_world, fname, "r")
+
+    mesh = Mesh()
+    subdomains = MeshFunction("size_t", mesh, 2)
+    surfaces = MeshFunction("size_t", mesh, 1)
+    hdf5file.read(mesh, '/mesh', False)
+    mesh.coordinates()[:] *= 1e6
+    hdf5file.read(subdomains, '/subdomains')
+    hdf5file.read(surfaces, '/surfaces')
+
+    x_min = 20.3; x_max = 20.4
+    y_min = 0.76; y_max = 0.77
+    z_min = 0.79; z_max = 0.81
+
+    # define one facet to 10 for getting membrane potential
+    for facet in facets(mesh):
+        x = [facet.midpoint().x(), facet.midpoint().y(), facet.midpoint().z()]
+
+        point_1 = x_min < x[0] < x_max and \
+                  y_min < x[1] < y_max and \
+                  z_min < x[2] < z_max
+
+        if point_1:
+            print("point", x[0], x[1], x[2])
+            surfaces[facet] = 10
+
+    surfacesfile = File('surfaces_plot.pvd')
+    surfacesfile << surfaces
+
+    # define function space of piecewise constants on interface gamma for solution to ODEs
+    Q = FunctionSpace(mesh, 'Discontinuous Lagrange Trace', 0)
+    phi_M = Function(Q)
+    E_Na = Function(Q)
+    E_K = Function(Q)
+    # interface normal
+    n_g = interface_normal(subdomains, mesh)
+
+    dS = Measure('dS', domain=mesh, subdomain_data=surfaces)
+    iface_size = assemble(Constant(1)*dS(10))
+
+    P1 = FiniteElement('DG', mesh.ufl_cell(), 1)
+    W = FunctionSpace(mesh, MixedElement(2*[P1]))
+    V = FunctionSpace(mesh, P1)
+
+    v = Function(W)
+    w_Na = Function(V)
+    w_phi = Function(V)
+
+    f_phi = Function(V)
+    f_Na = Function(V)
+    f_K = Function(V)
+
+    phi_M_s = []
+    E_Na_s = []
+    E_K_s = []
+
+    z_Na = 1; z_K = 1; temperature = 300; F = 96485; R = 8.314
+
+    for n in range(1, int(T/dt)):
+
+            # potential
+            hdf5file.read(w_phi, "/potential/vector_" + str(n))
+            assign(f_phi, w_phi)
+
+            # K concentrations
+            hdf5file.read(v, "/concentrations/vector_" + str(n))
+            assign(f_K, v.sub(0))
+            E = R * temperature / (F * z_K) * ln(plus(f_K, n_g) / minus(f_K, n_g))
+            assign(E_K, pcws_constant_project(E, Q))
+            E_K_ = 1.0e3*assemble(1.0/iface_size*avg(E_K)*dS(10))
+            E_K_s.append(E_K_)
+
+            # Na concentrations
+            hdf5file.read(w_Na, "/elim_concentration/vector_" + str(n))
+            assign(f_Na, w_Na)
+            E = R * temperature / (F * z_Na) * ln(plus(f_Na, n_g) / minus(f_Na, n_g))
+            assign(E_Na, pcws_constant_project(E, Q))
+            E_Na_ = 1.0e3*assemble(1.0/iface_size*avg(E_Na)*dS(10))
+            E_Na_s.append(E_Na_)
+
+            # update membrane potential
+            phi_M_step = JUMP(f_phi, n_g)
+            assign(phi_M, pcws_constant_project(phi_M_step, Q))
+            phi_M_s.append(1.0e3*assemble(1.0/iface_size*avg(phi_M)*dS(10)))
+
+    return phi_M_s, E_Na_s, E_K_s
+
+
+def plot_3D_concentration(res, T, dt):
+
+    temperature = 300 # temperature (K)
+    F = 96485         # Faraday's constant (C/mol)
+    R = 8.314         # Gas constant (J/(K*mol))
+
+    time = 1.0e3*np.arange(0, T-dt, dt)
+
+    # at membrane of axon A (gamma)
+    x_M_A = 25.6; y_M_A = 0.333; z_M_A = 0.8
+    # 0.05 um above axon A (ECS)
+    x_e_A = 25; y_e_A = 0.9; z_e_A = 0.9
+    #x_e_A = 50; y_e_A = 1.0; z_e_A = 0.7
+    # mid point inside axon A (ICS)
+    x_i_A = 25; y_i_A = 0.3; z_i_A = 0.6
+
+    #################################################################
+    # get data axon A is stimulated
+    fname = 'results/data/3D/results.h5'
+
+    # trace concentrations
+    phi_M, E_Na, E_K = get_time_series_membrane(dt, T, fname, x_M_A, y_M_A, z_M_A)
+
+    # bulk concentrations
+    Na_e, K_e, Cl_e, _ = get_time_series(dt, T, fname, x_e_A, y_e_A, z_e_A)
+    Na_i, K_i, Cl_i, _ = get_time_series(dt, T, fname, x_i_A, y_i_A, z_i_A)
+
+    #################################################################
+    # get data axons BC are stimulated
+
+    # Concentration plots
+    fig = plt.figure(figsize=(12*0.9,12*0.9))
+    ax = plt.gca()
+
+    print("-------------------------")
+    print("Na_e", Na_e[-1])
+    print("K_e", K_e[-1])
+    print("Cl_e", Cl_e[-1])
+
+    print("Na_i", Na_i[-1])
+    print("K_i", K_i[-1])
+    print("Cl_i", Cl_i[-1])
+
+    print("phi_M", phi_M[-1])
+    print("-------------------------")
+
+    ax1 = fig.add_subplot(3,3,1)
+    plt.title(r'Na$^+$ concentration (ECS)')
+    plt.ylabel(r'[Na]$_e$ (mM)')
+    plt.plot(Na_e, linewidth=3, color='b')
+
+    ax3 = fig.add_subplot(3,3,2)
+    plt.title(r'K$^+$ concentration (ECS)')
+    plt.ylabel(r'[K]$_e$ (mM)')
+    plt.plot(K_e, linewidth=3, color='b')
+
+    ax3 = fig.add_subplot(3,3,3)
+    plt.title(r'Cl$^-$ concentration (ECS)')
+    plt.ylabel(r'[Cl]$_e$ (mM)')
+    plt.plot(Cl_e, linewidth=3, color='b')
+
+    ax2 = fig.add_subplot(3,3,4)
+    plt.title(r'Na$^+$ concentration (ICS)')
+    plt.ylabel(r'[Na]$_i$ (mM)')
+    plt.plot(Na_i,linewidth=3, color='r')
+
+    ax2 = fig.add_subplot(3,3,5)
+    plt.title(r'K$^+$ concentration (ICS)')
+    plt.ylabel(r'[K]$_i$ (mM)')
+    plt.plot(K_i,linewidth=3, color='r')
+
+    ax2 = fig.add_subplot(3,3,6)
+    plt.title(r'Cl$^-$ concentration (ICS)')
+    plt.ylabel(r'[Cl]$_i$ (mM)')
+    plt.plot(Cl_i,linewidth=3, color='r')
+
+    ax5 = fig.add_subplot(3,3,7)
+    plt.title(r'Membrane potential')
+    plt.ylabel(r'$\phi_M$ (mV)')
+    plt.xlabel(r'time (ms)')
+    plt.plot(phi_M, linewidth=3)
+
+    ax6 = fig.add_subplot(3,3,8)
+    plt.title(r'Na$^+$ reversal potential')
+    plt.ylabel(r'E$_Na$ (mV)')
+    plt.xlabel(r'time (ms)')
+    plt.plot(E_K, linewidth=3)
+    plt.plot(E_Na, linewidth=3)
+
+    # make pretty
+    ax.axis('off')
+    plt.tight_layout()
+
+    # save figure to file
+    plt.savefig('results/figures/pot_con_3D.svg', format='svg')
+
+    f_phi_M = open('results/data/3D/solver/phi_M_3D.txt', "w")
+    for p in phi_M:
+        f_phi_M.write("%.10f \n" % p*1000)
+    f_phi_M.close()
+
+    return
+
+def get_plottable_gamma_function(fname, n):
+    """ Return list of values in given point (x, y, z) over time """
+
+    # read data file
+    hdf5file = HDF5File(MPI.comm_world, fname, "r")
+
+    # create mesh
+    mesh = Mesh()
+    subdomains = MeshFunction("size_t", mesh, mesh.topology().dim())
+    surfaces = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
+    hdf5file.read(mesh, "/mesh", False)
+    mesh.coordinates()[:] *= 1e6
+    hdf5file.read(subdomains, "/subdomains")
+    hdf5file.read(surfaces, "/surfaces")
+
+    mesh = mesh
+
+    # create function spaces
+    P1 = FiniteElement("DG", mesh.ufl_cell(), 1)
+    V = FunctionSpace(mesh, P1)
+
+    # create functions
+    f_phi = Function(V)
+    w = Function(V)
+
+    # potential
+    hdf5file.read(w, "/potential/vector_" + str(n))
+    assign(f_phi, w)
+
+    # membrane potential
+    Q = FunctionSpace(mesh, 'Discontinuous Lagrange Trace', 0)
+    dS = Measure('dS', domain=mesh, subdomain_data=surfaces)
+    n_g = interface_normal(subdomains, mesh)
+    phi_M = Function(Q)
+
+    # update membrane potential
+    phi_M_step = JUMP(f_phi, n_g)
+    assign(phi_M, pcws_constant_project(phi_M_step, Q))
+
+    phi_M_s = []
+
+    x_min = 1; x_max = 32
+    y_min = 0.76; y_max = 0.77
+    z_min = 0.79; z_max = 0.81
+
+    # define one facet to 10 for getting membrane potential
+    for facet in facets(mesh):
+        x = [facet.midpoint().x(), facet.midpoint().y(), facet.midpoint().z()]
+        line = (x_min <= x[0] <= x_max and y_min <= x[1] <= y_max and  z_min <= x[2] <= z_max)
+        if (surfaces[facet] == 1 and line):
+            print(x[0], x[1], x[2])
+            surfaces[facet] = 10
+
+            iface_size = assemble(Constant(1)*dS(10))
+            phi_M_s.append(1.0e3*assemble(1.0/iface_size*avg(phi_M)*dS(10)))
+
+        if (surfaces[facet] == 1 and line):
+            print(x[0], x[1], x[2])
+
+    return phi_M_s
+
+def plot_membrane_space(fname, n):
+
+    plt.figure(figsize=(4.6, 3.3))
+    plt.xlabel(r"x-position ($\mu$m)")
+    plt.ylabel("$\phi_M$ (mV)")
+    plt.ylim([-80, -60])
+    #plt.ylim([-100e-3, -70e-3])
+    #plt.yticks([-80, -60, -40, -20, 0, 20, 40])
+    #plt.ylim([-90, 50])
+
+    colors_n = [c_1, c_2, c_3]
+
+    header = "x r0 r1 r2"
+    phi_dat = []
+
+    for idx, n in enumerate([1, 50, 99]):
+        phi_M_n = get_plottable_gamma_function(fname, n)
+        n_new = n/10. # convert n to ms
+        x = np.linspace(10, 190, len(phi_M_n))
+        plt.plot(x, phi_M_n, linewidth=2.5, color=colors_n[idx], label=r"%d ms" % n_new)
+        phi_dat.append(phi_M_n)
+
+    plt.legend()
+    #plt.tight_layout()
+    # save figure
+    plt.savefig("results/figures/phi_M_space.png")
+    plt.close()
+
+    phi_dat.insert(0, x)
+    a_phi = np.asarray(phi_dat)
+    np.savetxt("results/figures/phi_space.dat", np.transpose(a_phi), delimiter=" ", header=header)
+
+    return
+
+
+# create directory for figures
+if not os.path.isdir('results/figures'):
+    os.mkdir('results/figures')
+
+# create figures
+res_3D = '0' # mesh resolution for 3D axon bundle
+#T = 1.0e-1
+T = 5.0e-2
+dt = 1.0e-4
+
+n = 99
+fname = 'results/data/3D/results.h5'
+
+#plot_3D_concentration(res_3D, T, dt)
+plot_membrane_space(fname, n)
diff --git a/elec-prop-test/make_figures_3D_emi.py b/elec-prop-test/make_figures_3D_emi.py
new file mode 100644
index 0000000..1e726c1
--- /dev/null
+++ b/elec-prop-test/make_figures_3D_emi.py
@@ -0,0 +1,496 @@
+import matplotlib as mpl
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+import sys
+
+from fenics import * 
+import string
+
+from knpemidg import pcws_constant_project
+from knpemidg import interface_normal, plus, minus
+
+JUMP = lambda f, n: minus(f, n) - plus(f, n)
+
+# set font & text parameters
+font = {'family' : 'serif',
+        'weight' : 'bold',
+        'size'   : 13}
+
+plt.rc('font', **font)
+plt.rc('text', usetex=True)
+mpl.rcParams['image.cmap'] = 'jet'
+
+path = 'results/data/'
+
+c_1 = '#%02x%02x%02x' % (255, 0, 255)  # pink
+c_2 = '#%02x%02x%02x' % (54, 92, 141)  # blue
+c_3 = '#%02x%02x%02x' % (255, 165, 0)  # orange
+phi_M_c = '#%02x%02x%02x' % (54, 92, 141)
+
+def get_time_series(dt, T, fname, x, y, z):
+    # read data file
+    hdf5file = HDF5File(MPI.comm_world, fname, "r")
+
+    mesh = Mesh()
+    subdomains = MeshFunction("size_t", mesh, 2)
+    surfaces = MeshFunction("size_t", mesh, 1)
+    hdf5file.read(mesh, '/mesh', False)
+    mesh.coordinates()[:] *= 1e4
+    hdf5file.read(subdomains, '/subdomains')
+    hdf5file.read(surfaces, '/surfaces')
+
+    P1 = FiniteElement('CG', mesh.ufl_cell(), 1)
+    W = FunctionSpace(mesh, MixedElement(2*[P1]))
+    V = FunctionSpace(mesh, P1)
+
+    u = Function(W)
+    v = Function(V)
+    w = Function(V)
+
+    f_Na = Function(V)
+    f_K = Function(V)
+    f_Cl = Function(V)
+    f_phi = Function(V)
+
+    Na = []
+    K = []
+    Cl = []
+    phi = []
+
+    for n in range(1, int(T/dt)):
+
+            # read file
+            hdf5file.read(u, "/concentrations/vector_" + str(n))
+
+            # K concentrations
+            assign(f_K, u.sub(0))
+            K.append(f_K(x, y, z))
+
+            # Cl concentrations
+            assign(f_Cl, u.sub(1))
+            Cl.append(f_Cl(x, y, z))
+
+            # Na concentrations
+            hdf5file.read(v, "/elim_concentration/vector_" + str(n))
+            assign(f_Na, v)
+            Na.append(f_Na(x, y, z))
+
+            # potential
+            hdf5file.read(w, "/potential/vector_" + str(n))
+            assign(f_phi, w)
+            phi.append(f_phi(x, y, z))
+
+    return Na, K, Cl, phi
+
+def get_time_series_membrane(dt, T, fname, x_, y_, z_):
+    # read data file
+    hdf5file = HDF5File(MPI.comm_world, fname, "r")
+
+    mesh = Mesh()
+    subdomains = MeshFunction("size_t", mesh, 2)
+    surfaces = MeshFunction("size_t", mesh, 1)
+    hdf5file.read(mesh, '/mesh', False)
+    mesh.coordinates()[:] *= 1e4
+    hdf5file.read(subdomains, '/subdomains')
+    hdf5file.read(surfaces, '/surfaces')
+
+    x_min = 100.3; x_max = 100.4
+    y_min = 0.76; y_max = 0.77
+    z_min = 0.79; z_max = 0.81
+
+    # define one facet to 10 for getting membrane potential
+    for facet in facets(mesh):
+        x = [facet.midpoint().x(), facet.midpoint().y(), facet.midpoint().z()]
+
+        point_1 = x_min < x[0] < x_max and \
+                  y_min < x[1] < y_max and \
+                  z_min < x[2] < z_max
+
+        if point_1:
+            surfaces[facet] = 10
+
+    surfacesfile = File('surfaces_plot.pvd')
+    surfacesfile << surfaces
+
+    # define function space of piecewise constants on interface gamma for solution to ODEs
+    Q = FunctionSpace(mesh, 'Discontinuous Lagrange Trace', 0)
+    phi_M = Function(Q)
+    E_Na = Function(Q)
+    E_K = Function(Q)
+
+    # interface normal
+    n_g = interface_normal(subdomains, mesh)
+
+    dS = Measure('dS', domain=mesh, subdomain_data=surfaces)
+    iface_size = assemble(Constant(1)*dS(10))
+
+    P1 = FiniteElement('DG', mesh.ufl_cell(), 1)
+    W = FunctionSpace(mesh, MixedElement(2*[P1]))
+    V = FunctionSpace(mesh, P1)
+
+    v = Function(W)
+    w_Na = Function(V)
+    w_phi = Function(V)
+
+    f_phi = Function(V)
+    f_Na = Function(V)
+    f_K = Function(V)
+
+    phi_M_s = []
+    E_Na_s = []
+    E_K_s = []
+
+    z_Na = 1; z_K = 1; temperature = 300; F = 96.485; R = 8.314
+
+    for n in range(1, int(T/dt)):
+
+            # potential
+            hdf5file.read(w_phi, "/potential/vector_" + str(n))
+            assign(f_phi, w_phi)
+
+            # K concentrations
+            hdf5file.read(v, "/concentrations/vector_" + str(n))
+            assign(f_K, v.sub(0))
+            E = R * temperature / (F * z_K) * ln(plus(f_K, n_g) / minus(f_K, n_g))
+            assign(E_K, pcws_constant_project(E, Q))
+            E_K_ = assemble(1.0/iface_size*avg(E_K)*dS(10))
+            E_K_s.append(E_K_)
+
+            # Na concentrations
+            hdf5file.read(w_Na, "/elim_concentration/vector_" + str(n))
+            assign(f_Na, w_Na)
+            E = R * temperature / (F * z_Na) * ln(plus(f_Na, n_g) / minus(f_Na, n_g))
+            assign(E_Na, pcws_constant_project(E, Q))
+            E_Na_ = assemble(1.0/iface_size*avg(E_Na)*dS(10))
+            E_Na_s.append(E_Na_)
+
+            # update membrane potential
+            phi_M_step = JUMP(f_phi, n_g)
+            assign(phi_M, pcws_constant_project(phi_M_step, Q))
+            phi_M_s.append(assemble(1.0/iface_size*avg(phi_M)*dS(10)))
+
+    return phi_M_s, E_Na_s, E_K_s
+
+
+def plot_3D_concentration(res, T, dt):
+
+    time = np.arange(0, T-dt, dt)
+
+    # at membrane of axon A (gamma)
+    x_M_A = 100; y_M_A = 0.333; z_M_A = 0.8
+    # 0.05 um above axon A (ECS)
+    x_e_A = 100; y_e_A = 0.9; z_e_A = 0.9
+    # mid point inside axon A (ICS)
+    x_i_A = 100; y_i_A = 0.3; z_i_A = 0.6
+
+    #################################################################
+    # get data axon A is stimulated
+    fname = 'results/data/3D_emi/results.h5'
+
+    # trace concentrations
+    phi_M, E_Na, E_K = get_time_series_membrane(dt, T, fname, x_M_A, y_M_A, z_M_A)
+
+    # bulk concentrations
+    Na_e, K_e, Cl_e, _ = get_time_series(dt, T, fname, x_e_A, y_e_A, z_e_A)
+    Na_i, K_i, Cl_i, _ = get_time_series(dt, T, fname, x_i_A, y_i_A, z_i_A)
+
+    #################################################################
+    # get data axons BC are stimulated
+
+    # Concentration plots
+    fig = plt.figure(figsize=(12*0.9,12*0.9))
+    ax = plt.gca()
+
+    print("-------------------------")
+    print("Na_e", Na_e[-1])
+    print("K_e", K_e[-1])
+    print("Cl_e", Cl_e[-1])
+
+    print("Na_i", Na_i[-1])
+    print("K_i", K_i[-1])
+    print("Cl_i", Cl_i[-1])
+
+    print("phi_M", phi_M[-1])
+    print("E_Na", E_Na[-1])
+    print("E_K", E_K[-1])
+    print("-------------------------")
+
+    # global kappa
+    Na_i_init = Na_i[-1]
+    Na_e_init = Na_e[-1]
+    K_i_init = K_i[-1]
+    K_e_init = K_e[-1]
+    Cl_i_init = Cl_i[-1]
+    Cl_e_init = Cl_e[-1]
+
+    # Physical parameters (PDEs)
+    C_M = 1.0                      # capacitance
+    temperature = 300e3            # temperature (m K)
+    R = 8.314e3                    # Gas Constant (m J/(K mol))
+    F = 96485e3                    # Faraday's constant (mC/ mol)
+
+    D_Na = 1.33e-8       # diffusion coefficients Na (cm/ms)
+    D_K = 1.96e-8        # diffusion coefficients K (cm/ms)
+    D_Cl = 2.03e-8       # diffusion coefficients Cl (cm/ms)
+
+    psi = F/(R * temperature)
+
+    kappa_i = F * psi * (float(D_Na) * float(Na_i_init) \
+                      + float(D_K) * float(K_i_init) \
+                      + float(D_Cl) * float(Cl_i_init))
+
+    kappa_e = F * psi * (float(D_Na) * float(Na_e_init) \
+                       + float(D_K) * float(K_e_init) \
+                       + float(D_Cl) * float(Cl_e_init))
+
+    E_K = float((R * temperature)/F)
+
+    print("kappa_i", kappa_i)
+    print("kappa_e", kappa_e)
+
+    print(kappa_e)
+    print(E_K)
+    print(float((F*F)/(R*temperature)))
+
+    ax1 = fig.add_subplot(3,3,1)
+    plt.title(r'Na$^+$ concentration (ECS)')
+    plt.ylabel(r'[Na]$_e$ (mM)')
+    plt.plot(Na_e, linewidth=3, color='b')
+
+    ax3 = fig.add_subplot(3,3,2)
+    plt.title(r'K$^+$ concentration (ECS)')
+    plt.ylabel(r'[K]$_e$ (mM)')
+    plt.plot(K_e, linewidth=3, color='b')
+
+    ax3 = fig.add_subplot(3,3,3)
+    plt.title(r'Cl$^-$ concentration (ECS)')
+    plt.ylabel(r'[Cl]$_e$ (mM)')
+    plt.plot(Cl_e, linewidth=3, color='b')
+
+    ax2 = fig.add_subplot(3,3,4)
+    plt.title(r'Na$^+$ concentration (ICS)')
+    plt.ylabel(r'[Na]$_i$ (mM)')
+    plt.plot(Na_i,linewidth=3, color='r')
+
+    ax2 = fig.add_subplot(3,3,5)
+    plt.title(r'K$^+$ concentration (ICS)')
+    plt.ylabel(r'[K]$_i$ (mM)')
+    plt.plot(K_i,linewidth=3, color='r')
+
+    ax2 = fig.add_subplot(3,3,6)
+    plt.title(r'Cl$^-$ concentration (ICS)')
+    plt.ylabel(r'[Cl]$_i$ (mM)')
+    plt.plot(Cl_i,linewidth=3, color='r')
+
+    ax5 = fig.add_subplot(3,3,7)
+    plt.title(r'Membrane potential')
+    plt.ylabel(r'$\phi_M$ (mV)')
+    plt.xlabel(r'time (ms)')
+    plt.plot(phi_M, linewidth=3)
+
+    ax6 = fig.add_subplot(3,3,8)
+    plt.title(r'Na$^+$ reversal potential')
+    plt.ylabel(r'E$_Na$ (mV)')
+    plt.xlabel(r'time (ms)')
+    plt.plot(E_K, linewidth=3)
+    plt.plot(E_Na, linewidth=3)
+
+    # make pretty
+    ax.axis('off')
+    plt.tight_layout()
+
+    # save figure to file
+    plt.savefig('results/figures/pot_con_3D_emi.svg', format='svg')
+
+    f_phi_M = open('results/data/3D_emi/solver/phi_M_3D.txt', "w")
+    for p in phi_M:
+        f_phi_M.write("%.10f \n" % p)
+    f_phi_M.close()
+
+    return
+
+def get_plottable_gamma_function(fname, n):
+    """ Return list of values in given point (x, y, z) over time """
+
+    # read data file
+    hdf5file = HDF5File(MPI.comm_world, fname, "r")
+
+    # create mesh
+    mesh = Mesh()
+    subdomains = MeshFunction("size_t", mesh, mesh.topology().dim())
+    surfaces = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
+    hdf5file.read(mesh, "/mesh", False)
+    mesh.coordinates()[:] *= 1e4
+    hdf5file.read(subdomains, "/subdomains")
+    hdf5file.read(surfaces, "/surfaces")
+
+    mesh = mesh
+
+    # create function spaces
+    P1 = FiniteElement("DG", mesh.ufl_cell(), 1)
+    V = FunctionSpace(mesh, P1)
+
+    # create functions
+    f_phi = Function(V)
+    w = Function(V)
+
+    # potential
+    hdf5file.read(w, "/potential/vector_" + str(n))
+    assign(f_phi, w)
+
+    # membrane potential
+    Q = FunctionSpace(mesh, 'Discontinuous Lagrange Trace', 0)
+    dS = Measure('dS', domain=mesh, subdomain_data=surfaces)
+    n_g = interface_normal(subdomains, mesh)
+    phi_M = Function(Q)
+
+    # update membrane potential
+    phi_M_step = JUMP(f_phi, n_g)
+    assign(phi_M, pcws_constant_project(phi_M_step, Q))
+
+    phi_M_s = []
+    x_s = []
+
+    x_min = 5.0; x_max = 205.0
+    y_min = 0.76; y_max = 0.77
+    z_min = 0.79; z_max = 0.81
+
+    # define one facet to 10 for getting membrane potential
+    for facet in facets(mesh):
+        x = [facet.midpoint().x(), facet.midpoint().y(), facet.midpoint().z()]
+        line = (x_min <= x[0] <= x_max and y_min <= x[1] <= y_max and  z_min <= x[2] <= z_max)
+        if (surfaces[facet] == 1 and line):
+            print(x[0], x[1], x[2])
+            surfaces[facet] = 10
+
+            iface_size = assemble(Constant(1)*dS(10))
+            phi_M_s.append(assemble(1.0/iface_size*avg(phi_M)*dS(10)))
+            x_s.append(x[0])
+
+        #if (surfaces[facet] == 1):
+            #print(x[0], x[1], x[2])
+
+    return phi_M_s, x_s
+
+def plot_membrane_space(fname):
+
+    plt.figure(figsize=(4.6, 3.3))
+    plt.xlabel(r"x-position ($\mu$m)")
+    plt.ylabel("$\phi_M$ (mV)")
+
+    header = "x r0 r1 r2"
+    phi_dat = []
+
+    for idx, n in enumerate([1, 5, 10, 15, 20, 25, 30, 35, 99]):
+        phi_M_n, x_s = get_plottable_gamma_function(fname, n)
+        n_new = n/10
+        x = np.linspace(5, 200, len(phi_M_n))
+        plt.plot(x_s, phi_M_n, linewidth=2.5, label=r"%d ms" % n_new)
+        phi_dat.append(phi_M_n)
+
+    plt.legend()
+    plt.tight_layout()
+    # save figure
+    plt.savefig("results/figures/phi_M_space_emi.svg")
+    plt.close()
+
+    #phi_dat.insert(0, x)
+    #a_phi = np.asarray(phi_dat)
+    #np.savetxt("results/figures/phi_space_emi.dat", np.transpose(a_phi), delimiter=" ", header=header)
+
+    return
+
+def get_velocity(fname, T, dt):
+    # read data file
+    hdf5file = HDF5File(MPI.comm_world, fname, "r")
+
+    mesh = Mesh()
+    subdomains = MeshFunction("size_t", mesh, 2)
+    surfaces = MeshFunction("size_t", mesh, 1)
+    hdf5file.read(mesh, '/mesh', False)
+    mesh.coordinates()[:] *= 1e4
+    hdf5file.read(subdomains, '/subdomains')
+    hdf5file.read(surfaces, '/surfaces')
+
+    # define one facet to 10 for getting membrane potential
+    for facet in facets(mesh):
+        x = [facet.midpoint().x(), facet.midpoint().y(), facet.midpoint().z()]
+        # mark top point
+        if surfaces[facet] == 1 and (20 < x[0] < 20.5) and (x[1] > 0.79) and (x[2] > 0.76):
+            print("point 1", x[0], x[1], x[2])
+            surfaces[facet] = 10
+        # mark bottom point
+        if surfaces[facet] == 1 and (120 < x[0] < 120.5) and (x[1] > 0.79) and (x[2] > 0.76):
+            print("point 2", x[0], x[1], x[2])
+            surfaces[facet] = 20
+
+    # define function space of piecewise constants on interface gamma for solution to ODEs
+    Q = FunctionSpace(mesh, 'Discontinuous Lagrange Trace', 0)
+    phi_M = Function(Q)
+
+    # interface normal
+    n_g = interface_normal(subdomains, mesh)
+
+    dS = Measure('dS', domain=mesh, subdomain_data=surfaces)
+    iface_size_10 = assemble(Constant(1)*dS(10))
+    iface_size_20 = assemble(Constant(1)*dS(20))
+
+    P1 = FiniteElement('DG', mesh.ufl_cell(), 1)
+    V = FunctionSpace(mesh, P1)
+
+    w_phi = Function(V)
+    f_phi = Function(V)
+
+    phi_M_s_1 = []
+    phi_M_s_2 = []
+
+    time_point_1 = 0
+    time_point_2 = 0
+
+    for n in range(1, int(T/dt)):
+
+            # potential
+            hdf5file.read(w_phi, "/potential/vector_" + str(n))
+            assign(f_phi, w_phi)
+
+            # update membrane potential
+            phi_M_step = JUMP(f_phi, n_g)
+            assign(phi_M, pcws_constant_project(phi_M_step, Q))
+
+            if 1.0e3*assemble(1.0/iface_size_10*avg(phi_M)*dS(10)) > 20 and time_point_1 == 0:
+                time_point_1 = n*0.1
+
+            if 1.0e3*assemble(1.0/iface_size_10*avg(phi_M)*dS(20)) > 20 and time_point_2 == 0:
+                time_point_2 = n*0.1
+
+            print("t1", time_point_1)
+            print("t2", time_point_2)
+
+            # if membrane potential has reach 0 in both points then break
+            if (time_point_1 > 0) and (time_point_2 > 0):
+                break
+
+    delta_t = (time_point_2 - time_point_1) # ms
+    delta_x = 100 # um
+
+    print("velocity (um/ms)", delta_x/delta_t)
+
+    return phi_M_s_1, phi_M_s_2
+
+# create directory for figures
+if not os.path.isdir('results/figures'):
+    os.mkdir('results/figures')
+
+# create figures
+res_3D = '0' # mesh resolution for 3D axon bundle
+T = 30
+dt = 0.1
+
+fname = 'results/data/3D_emi/results.h5'
+
+plot_3D_concentration(res_3D, T, dt)
+plot_membrane_space(fname)
+get_velocity(fname, T, dt)
diff --git a/elec-prop-test/make_mesh_3D.py b/elec-prop-test/make_mesh_3D.py
new file mode 100755
index 0000000..2f970af
--- /dev/null
+++ b/elec-prop-test/make_mesh_3D.py
@@ -0,0 +1,117 @@
+#!/usr/bin/env python3
+
+"""
+This script generates a 3D mesh representing 4 axons.
+"""
+
+from dolfin import *
+import sys
+
+class Boundary(SubDomain):
+    # define exterior boundary
+    def inside(self, x, on_boundary):
+        return on_boundary
+
+def add_axon(mesh, subdomains, surfaces, a, b):
+    # define interior domain
+    in_interior = """ (x[0] >= %g && x[0] <= %g &&
+                       x[1] >= %g && x[1] <= %g &&
+                       x[2] >= %g && x[2] <= %g) """ % (
+        a[0],
+        b[0],
+        a[1],
+        b[1],
+        a[2],
+        b[2],
+    )
+
+    interior = CompiledSubDomain(in_interior)
+
+    # mark interior and exterior domain
+    for cell in cells(mesh):
+        x = cell.midpoint().array()
+        if (int(interior.inside(x, False))):
+            subdomains[cell] = 1
+    assert sum(1 for _ in SubsetIterator(subdomains, 1)) > 0
+
+    for facet in facets(mesh):
+        x = [facet.midpoint().x(), facet.midpoint().y(), facet.midpoint().z()]
+        side_1 = near(x[0], a[0]) and a[1] <= x[1] <= b[1] and a[2] <= x[2] <= b[2]
+        side_2 = near(x[0], b[0]) and a[1] <= x[1] <= b[1] and a[2] <= x[2] <= b[2]
+        side_3 = near(x[1], a[1]) and a[0] <= x[0] <= b[0] and a[2] <= x[2] <= b[2]
+        side_4 = near(x[1], b[1]) and a[0] <= x[0] <= b[0] and a[2] <= x[2] <= b[2]
+        side_5 = near(x[2], a[2]) and a[0] <= x[0] <= b[0] and a[1] <= x[1] <= b[1]
+        side_6 = near(x[2], b[2]) and a[0] <= x[0] <= b[0] and a[1] <= x[1] <= b[1]
+        surfaces[facet] += side_1 or side_2 or side_3 or side_4 or side_5 or side_6
+
+    return
+
+import argparse
+from pathlib import Path
+
+
+class CustomParser(
+    argparse.ArgumentDefaultsHelpFormatter, argparse.RawDescriptionHelpFormatter
+): ...
+
+def main(argv=None):
+    parser = argparse.ArgumentParser(formatter_class=CustomParser)
+    parser.add_argument(
+        "-r",
+        "--resolution",
+        dest="resolution_factor",
+        default=0,
+        type=int,
+        help="Mesh resolution factor",
+    )
+    parser.add_argument(
+        "-d",
+        "--directory",
+        dest="mesh_dir",
+        type=Path,
+        default=Path("meshes/MMS"),
+        help="Directory to save the mesh",
+    )
+    args = parser.parse_args(argv)
+    resolution_factor = args.resolution_factor
+    out_dir = args.mesh_dir
+
+    l = 2
+    nx = l * 105 * 2 ** resolution_factor
+    ny = 10 * 2 ** resolution_factor
+    nz = 10 * 2 ** resolution_factor
+
+    # box mesh
+    mesh = BoxMesh(Point(0, 0.0, 0.0), Point(l * 105, 1.0, 1.0), nx, ny, nz)
+    subdomains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
+    surfaces = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
+
+    a = Point(5, 0.2, 0.2)
+    b = Point(l * 105 - 5, 0.8, 0.8)
+    add_axon(mesh, subdomains, surfaces, a, b)
+
+    # mark exterior boundary
+    Boundary().mark(surfaces, 5)
+
+    # convert mesh from um to cm
+    mesh.coordinates()[:, :] *= 1e-4
+
+    # save .xml files
+    mesh_file = File(str((out_dir / f"mesh_{resolution_factor}.xml").absolute()))
+    mesh_file << mesh
+
+    subdomains_file = File(str((out_dir / f"subdomains_{resolution_factor}.xml").absolute()))
+    subdomains_file << subdomains
+
+    surfaces_file = File(str((out_dir / f"surfaces_{resolution_factor}.xml").absolute()))
+    surfaces_file << surfaces
+
+    # save .pvd files
+    mesh_file = File(str((out_dir / f"mesh_{resolution_factor}.pvd").absolute()))
+    mesh_file << mesh
+
+    subdomains_file = File(str((out_dir / f"subdomains_{resolution_factor}.pvd").absolute()))
+    subdomains_file << subdomains
+
+    surfaces_file = File(str((out_dir / f"surfaces_{resolution_factor}.pvd").absolute()))
+    surfaces_file << surfaces
diff --git a/elec-prop-test/mm_hh.py b/elec-prop-test/mm_hh.py
new file mode 100644
index 0000000..3e878f0
--- /dev/null
+++ b/elec-prop-test/mm_hh.py
@@ -0,0 +1,168 @@
+# Gotran generated code for the "hodgkin_huxley_squid_axon_model_1952_original" model
+
+import numpy as np
+import math
+
+def init_state_values(**values):
+    """
+    Initialize state values
+    """
+    # Init values
+    n_init = 0.27622914792          # gating variable n
+    m_init = 0.0379183462722        # gating variable m
+    h_init = 0.688489218108         # gating variable h
+    phi_M_init = -68.31807655880984 # membrane potential (mV)
+
+    init_values = np.array([m_init, h_init, n_init, phi_M_init], dtype=np.float_)
+
+    # State indices and limit checker
+    state_ind = dict([("m", 0), ("h", 1), ("n", 2), ("V", 3)])
+
+    for state_name, value in values.items():
+        if state_name not in state_ind:
+            raise ValueError("{0} is not a state.".format(state_name))
+        ind = state_ind[state_name]
+
+        # Assign value
+        init_values[ind] = value
+
+    return init_values
+
+def init_parameter_values(**values):
+    """
+    Initialize parameter values
+    """
+
+    # Membrane parameters
+    g_Na_bar = 120           # Na max conductivity (mS/cm**2)
+    g_K_bar = 36             # K max conductivity (mS/cm**2)
+    g_leak_Na = 0.06         # Na leak conductivity (mS/cm**2)
+    g_leak_K  = 0.19         # K leak conductivity (mS/cm**2)
+
+    # Set initial parameter values
+    init_values = np.array([g_Na_bar, g_K_bar, \
+                            g_leak_Na, g_leak_K, \
+                            0, 0, 0, 0, \
+                            0, 0, 0, 0, 0], dtype=np.float_)
+
+    # Parameter indices and limit checker
+    param_ind = dict([("g_Na_bar", 0), ("g_K_bar", 1), \
+                      ("g_leak_Na", 2), ("g_leak_K", 3), \
+                      ("E_Na", 4), ("E_K", 5), \
+                      ("Cm", 6), ("stim_amplitude", 7), \
+                      ("I_ch_Na", 8), ("I_ch_K", 9), ("I_ch_Cl", 10), \
+                      ("K_e", 11), ("Na_i", 12)])
+
+    for param_name, value in values.items():
+        if param_name not in param_ind:
+            raise ValueError("{0} is not a parameter.".format(param_name))
+        ind = param_ind[param_name]
+
+        # Assign value
+        init_values[ind] = value
+
+    return init_values
+
+def state_indices(*states):
+    """
+    State indices
+    """
+    state_inds = dict([("m", 0), ("h", 1), ("n", 2), ("V", 3)])
+
+    indices = []
+    for state in states:
+        if state not in state_inds:
+            raise ValueError("Unknown state: '{0}'".format(state))
+        indices.append(state_inds[state])
+    if len(indices)>1:
+        return indices
+    else:
+        return indices[0]
+
+def parameter_indices(*params):
+    """
+    Parameter indices
+    """
+    param_inds = dict([("g_Na_bar", 0), ("g_K_bar", 1), \
+                      ("g_leak_Na", 2), ("g_leak_K", 3), \
+                      ("E_Na", 4), ("E_K", 5), \
+                      ("Cm", 6), \
+                      ("stim_amplitude", 7), \
+                      ("I_ch_Na", 8), \
+                      ("I_ch_K", 9), \
+                      ("I_ch_Cl", 10), \
+                      ("K_e", 11), \
+                      ("Na_i", 12)])
+
+    indices = []
+    for param in params:
+        if param not in param_inds:
+            raise ValueError("Unknown param: '{0}'".format(param))
+        indices.append(param_inds[param])
+    if len(indices)>1:
+        return indices
+    else:
+        return indices[0]
+
+from numbalsoda import lsoda_sig
+from numba import njit, cfunc, jit
+import numpy as np
+import timeit
+import math
+
+@cfunc(lsoda_sig, nopython=True) 
+def rhs_numba(t, states, values, parameters):
+    """
+    Compute the right hand side of the\
+        hodgkin_huxley_squid_axon_model_1952_original ODE
+    """
+
+    # Assign states
+    #assert(len(states)) == 4
+
+    # Assign parameters
+    #assert(len(parameters)) == 11
+
+    # # Init return args
+    # if values is None:
+    #     values = np.zeros((4,), dtype=np.float_)
+    # else:
+    #     assert isinstance(values, np.ndarray) and values.shape == (4,)
+
+    alpha_m = 0.1 * (states[3] + 40.0)/(1.0 - math.exp(-(states[3] + 40.0) / 10.0))
+    beta_m = 4.0 * math.exp(-(states[3] + 65.0) / 18.0)
+
+    alpha_h = 0.07 * math.exp(-(states[3] + 65.0) / 20.0)
+    beta_h = 1.0 / (1.0 + math.exp(-(states[3] + 35.0) / 10.0))
+
+    alpha_n = 0.01 * (states[3] + 55.0)/(1.0 - math.exp(-(states[3] + 55.0) / 10.0))
+    beta_n = 0.125 * math.exp(-(states[3] + 65) / 80.0)
+
+    # Expressions for the m gate component
+    values[0] = (1 - states[0])*alpha_m - states[0]*beta_m
+
+    # Expressions for the h gate component
+    values[1] = (1 - states[1])*alpha_h - states[1]*beta_h
+
+    # Expressions for the n gate component
+    values[2] = (1 - states[2])*alpha_n - states[2]*beta_n
+
+    # Expressions for the Membrane component
+    i_Stim = parameters[7] * np.exp(-np.mod(t, 20.0)/2.0)
+
+    # Expressions for the Sodium channel component
+    i_Na = (parameters[2] + parameters[0]*states[1]*math.pow(states[0], 3) + i_Stim) * \
+           (states[3] - parameters[4])
+
+    # Expressions for the Potassium channel component
+    i_K = (parameters[3] + parameters[1]*math.pow(states[2], 4)) * \
+          (states[3] - parameters[5])
+
+    # set I_ch_Na
+    parameters[8] = i_Na
+    # set I_ch_K
+    parameters[9] = i_K
+    # set I_ch_Cl
+    parameters[10] = 0.0
+
+    values[3] = (- i_K - i_Na)/parameters[6]
diff --git a/elec-prop-test/mm_leak.py b/elec-prop-test/mm_leak.py
new file mode 100644
index 0000000..7e8a6ff
--- /dev/null
+++ b/elec-prop-test/mm_leak.py
@@ -0,0 +1,153 @@
+# Gotran generated code for the "hodgkin_huxley_squid_axon_model_1952_original" model
+
+import numpy as np
+import math
+
+def init_state_values(**values):
+    """
+    Initialize state values
+    """
+    # Init values
+    #phi_M_init = -0.06831807655880984
+    phi_M_init = -68.31807655880984
+
+    init_values = np.array([phi_M_init], dtype=np.float_)
+
+    # State indices and limit checker
+    state_ind = dict([("V", 0)])
+
+    for state_name, value in values.items():
+        if state_name not in state_ind:
+            raise ValueError("{0} is not a state.".format(state_name))
+        ind = state_ind[state_name]
+
+        # Assign value
+        init_values[ind] = value
+
+    return init_values
+
+def init_parameter_values(**values):
+    """
+    Initialize parameter values
+    """
+
+    # Membrane parameters
+    #g_leak_Na = 1.0     # Na leak conductivity (S/m**2)
+    #g_leak_K  = 4.0     # K leak conductivity (S/m**2)
+    #m_K = 4.0           # threshold ECS K (mol/m^3)
+    #m_Na = 12.0         # threshold ICS Na (mol/m^3)
+    #I_max = 1.3         # max pump strength (A/m^2)
+
+    g_leak_Na = 0.1     # Na leak conductivity (mS/cm**2)
+    g_leak_K  = 0.4     # K leak conductivity (mS/cm**2)
+
+    m_K = 4.0            # threshold ECS K (mol/m^3)
+    m_Na = 12.0          # threshold ICS Na (mol/m^3)
+    I_max = 1.3e-3       # max pump strength (muA/m^2)
+
+    # Set initial parameter values
+    init_values = np.array([g_leak_Na, g_leak_K, \
+                            0, 0, 0, 0, 0, 0, 0, \
+                            m_K, m_Na, I_max, \
+                            0, 0], dtype=np.float_)
+
+    # Parameter indices and limit checker
+    param_ind = dict([("g_leak_Na", 0), ("g_leak_K", 1), \
+                      ("E_Na", 2), ("E_K", 3), \
+                      ("Cm", 4), ("stim_amplitude", 5), \
+                      ("I_ch_Na", 6), ("I_ch_K", 7), ("I_ch_Cl", 8), \
+                      ("m_K", 9), ("m_Na", 10), ("I_max", 11), \
+                      ("K_e", 12), ("Na_i", 13)])
+
+    for param_name, value in values.items():
+        if param_name not in param_ind:
+            raise ValueError("{0} is not a parameter.".format(param_name))
+        ind = param_ind[param_name]
+
+        # Assign value
+        init_values[ind] = value
+
+    return init_values
+
+def state_indices(*states):
+    """
+    State indices
+    """
+    state_inds = dict([("V", 0)])
+
+    indices = []
+    for state in states:
+        if state not in state_inds:
+            raise ValueError("Unknown state: '{0}'".format(state))
+        indices.append(state_inds[state])
+    if len(indices)>1:
+        return indices
+    else:
+        return indices[0]
+
+def parameter_indices(*params):
+    """
+    Parameter indices
+    """
+    param_inds = dict([("g_leak_Na", 0), ("g_leak_K", 1), \
+                       ("E_Na", 2), ("E_K", 3), \
+                       ("Cm", 4), ("stim_amplitude", 5), \
+                       ("I_ch_Na", 6), ("I_ch_K", 7), ("I_ch_Cl", 8), \
+                       ("m_K", 9), ("m_Na", 10), ("I_max", 11), \
+                       ("K_e", 12), ("Na_i", 13)])
+
+    indices = []
+    for param in params:
+        if param not in param_inds:
+            raise ValueError("Unknown param: '{0}'".format(param))
+        indices.append(param_inds[param])
+    if len(indices)>1:
+        return indices
+    else:
+        return indices[0]
+
+from numbalsoda import lsoda_sig
+from numba import njit, cfunc, jit
+import numpy as np
+import timeit
+import math
+
+@cfunc(lsoda_sig, nopython=True) 
+def rhs_numba(t, states, values, parameters):
+    """
+    Compute the right hand side of the\
+        hodgkin_huxley_squid_axon_model_1952_original ODE
+    """
+
+    # Assign states
+    #assert(len(states)) == 4
+
+    # Assign parameters
+    #assert(len(parameters)) == 10
+
+    # Init return args
+    #if values is None:
+         #values = np.zeros((4,), dtype=np.float_)
+    #else:
+         #assert isinstance(values, np.ndarray) and values.shape == (4,)
+
+    # Expressions for stimuli
+    i_Stim = parameters[5]
+
+    I_pump = parameters[11] / ((1 + parameters[9] / parameters[12]) ** 2 \
+           * (1 + parameters[10] / parameters[13]) ** 3)
+
+    # Expressions for the Sodium channel component
+    i_Na = (parameters[0] + i_Stim) * (states[0] - parameters[2]) + 3 * I_pump
+
+    # Expressions for the Potassium channel component
+    i_K = parameters[1] * (states[0] - parameters[3]) - 2 * I_pump
+
+    # set I_ch_Na
+    parameters[6] = i_Na
+    # set I_ch_K
+    parameters[7] = i_K
+    # set I_ch_Cl
+    parameters[8] = 0.0
+
+    values[0] = (- i_K - i_Na)/parameters[4]
diff --git a/elec-prop-test/run_3D.py b/elec-prop-test/run_3D.py
new file mode 100644
index 0000000..3841bea
--- /dev/null
+++ b/elec-prop-test/run_3D.py
@@ -0,0 +1,155 @@
+#!/usr/bin/python3
+
+import os
+import sys
+import time
+
+from dolfin import *
+import numpy as np
+
+from collections import namedtuple
+
+from knpemidg import Solver
+import mm_leak as mm_leak
+
+# Define colors for printing
+class bcolors:
+    HEADER = '\033[95m'
+    OKBLUE = '\033[94m'
+    OKCYAN = '\033[96m'
+    OKGREEN = '\033[92m'
+    WARNING = '\033[93m'
+    FAIL = '\033[91m'
+    ENDC = '\033[0m'
+    BOLD = '\033[1m'
+    UNDERLINE = '\033[4m'
+
+if __name__ == "__main__":
+
+    # Resolution factor of mesh
+    #for resolution in [0, 1]:
+    for resolution in [0]:
+
+        # Time variables PDEs
+        dt = 1.0e-4                      # global time step (s)
+        Tstop = 5.0e-2                  # global end time (s)
+        t = Constant(0.0)                # time constant
+
+        # Time variables ODEs
+        n_steps_ODE = 25                 # number of ODE steps
+
+        # Physical parameters (PDEs)
+        C_M = 0.02                       # capacitance
+        temperature = 300                # temperature (K)
+        F = 96485                        # Faraday's constant (C/mol)
+        R = 8.314                        # Gas Constant (J/(K*mol))
+        D_Na = Constant(1.33e-9)         # diffusion coefficients Na (m/s)
+        D_K = Constant(1.96e-9)          # diffusion coefficients K (m/s)
+        D_Cl = Constant(2.03e-9)         # diffusion coefficients Cl (m/s)
+        psi = F / (R * temperature)      # shorthand
+        C_phi = C_M / dt                 # shorthand
+
+        Na_e_init = Constant(99.99996207620596)
+        K_e_init = Constant(4.000072443420079)
+        Cl_e_init = Constant(104.00003451962611)
+        Na_i_init = Constant(12.000067694452486)
+        K_i_init = Constant(124.9998827540323)
+        Cl_i_init = Constant(136.9999504484848)
+        phi_M_init = Constant(-0.06831807655880984)
+
+        # Set physical parameters
+        params = namedtuple('params', ('dt', 'n_steps_ODE', 'F', 'psi', \
+                    'phi_M_init', 'C_phi', 'C_M', 'R', 'temperature'))(dt, \
+                    n_steps_ODE, F, psi, phi_M_init, C_phi, C_M, R, temperature)
+
+        # diffusion coefficients for each sub-domain
+        D_Na_sub = {1:D_Na, 0:D_Na}
+        D_K_sub = {1:D_K, 0:D_K}
+        D_Cl_sub = {1:D_Cl, 0:D_Cl}
+
+        # initial concentrations for each sub-domain
+        Na_init_sub = {1:Na_i_init, 0:Na_e_init}
+        K_init_sub = {1:K_i_init, 0:K_e_init}
+        Cl_init_sub = {1:Cl_i_init, 0:Cl_e_init}
+
+        # Create ions (channel conductivity is set below in the membrane model)
+        Na = {'c_init_sub':Na_init_sub, 'bdry': Constant((0, 0)),
+              'z':1.0, 'name':'Na', 'D_sub':D_Na_sub}
+        K = {'c_init_sub':K_init_sub, 'bdry': Constant((0, 0)),
+             'z':1.0, 'name':'K', 'D_sub':D_K_sub}
+        Cl = {'c_init_sub':Cl_init_sub, 'bdry': Constant((0, 0)),
+              'z':-1.0, 'name':'Cl', 'D_sub':D_Cl_sub}
+
+        # Create ion list. NB! The last ion in list will be eliminated, and
+        # should be the ion with the smallest diffusion coefficient D_k
+        ion_list = [K, Cl, Na]
+
+        # Membrane parameters
+        g_syn_bar = 10                   # synaptic conductivity (S/m**2)
+
+        # Set stimulus ODE
+        stimulus = {'stim_amplitude': g_syn_bar}
+        stimulus_locator = lambda x: (x[0] < 20.0e-4)
+
+        stim_params = namedtuple('membrane_params', ('g_syn_bar', \
+                                 'stimulus', 'stimulus_locator'))(g_syn_bar, \
+                                  stimulus, stimulus_locator)
+
+        # Get mesh, subdomains, surfaces paths
+        mesh_prefix = "meshes/3D/"
+        mesh_path = mesh_prefix + "mesh_" + str(resolution) + ".xml"
+        subdomains_path = mesh_prefix + "subdomains_" + str(resolution) + ".xml"
+        surfaces_path = mesh_prefix + "surfaces_" + str(resolution) + ".xml"
+
+        # Generate mesh if it does not exist
+        if not os.path.isfile(mesh_path):
+            from make_mesh_3D import main
+            main(["-r", str(resolution), "-d", mesh_prefix])
+
+        mesh = Mesh(mesh_path)
+        subdomains = MeshFunction('size_t', mesh, subdomains_path)
+        surfaces = MeshFunction('size_t', mesh, surfaces_path)
+
+        # Set solver parameters EMI (True is direct, and False is iterate) 
+        direct_emi = False
+        rtol_emi = 1E-5
+        atol_emi = 1E-40
+        threshold_emi = 0.9
+
+        # Set solver parameters KNP (True is direct, and False is iterate) 
+        direct_knp = False
+        rtol_knp = 1E-7
+        atol_knp = 2E-40
+        threshold_knp = 0.75
+
+        # Set parameters
+        solver_params = namedtuple('solver_params', ('direct_emi',
+                                   'direct_knp', 'resolution',
+                                   'rtol_emi', 'rtol_knp',
+                                   'atol_emi', 'atol_knp',
+                                   'threshold_emi', 'threshold_knp'
+                                   ))(direct_emi, direct_knp, resolution, \
+                                      rtol_emi, rtol_knp, atol_emi, atol_knp, \
+                                      threshold_emi, threshold_knp)
+
+        # Set solver parameters (True is direct, and False is iterate)
+        #direct_emi = False
+        #direct_knp = False
+
+        # Set solver parameters
+        #solver_params = namedtuple('solver_params', ('direct_emi',
+            #'direct_knp', 'resolution'))(direct_emi, direct_knp, resolution)
+
+        # File for results
+        fname = "results/data/3D/"
+
+        # Dictionary with membrane models (key is facet tag, value is ode model)
+        ode_models = {1: mm_leak}
+
+        # Solve system
+        S = Solver(params, ion_list)                    # create solver
+        S.setup_domain(mesh, subdomains, surfaces)      # setup meshes
+        S.setup_parameters()                            # setup physical parameters
+        S.setup_FEM_spaces()                            # setup function spaces and numerical parameters
+        S.setup_membrane_model(stim_params, ode_models) # setup membrane model(s)
+        S.solve_system_active(Tstop, t, solver_params, filename=fname) # solve
diff --git a/elec-prop-test/run_3D_emi.py b/elec-prop-test/run_3D_emi.py
new file mode 100644
index 0000000..b67117e
--- /dev/null
+++ b/elec-prop-test/run_3D_emi.py
@@ -0,0 +1,171 @@
+#!/usr/bin/python3
+
+import os
+import sys
+import time
+
+from dolfin import *
+import numpy as np
+
+from collections import namedtuple
+
+from knpemidg import SolverEMI
+import mm_leak as mm_leak
+import mm_hh as mm_hh
+
+# Define colors for printing
+class bcolors:
+    HEADER = '\033[95m'
+    OKBLUE = '\033[94m'
+    OKCYAN = '\033[96m'
+    OKGREEN = '\033[92m'
+    WARNING = '\033[93m'
+    FAIL = '\033[91m'
+    ENDC = '\033[0m'
+    BOLD = '\033[1m'
+    UNDERLINE = '\033[4m'
+
+if __name__ == "__main__":
+
+    # Resolution factor of mesh
+    #for resolution in [0, 1]:
+    for resolution in [0]:
+
+        # Time variables PDEs
+        dt = 0.1                       # global time step (ms)
+        Tstop = 30                     # global end time (ms)
+        t = Constant(0.0)              # time constant
+
+        # Time variables ODEs
+        n_steps_ODE = 20               # number of ODE steps
+
+        # Physical parameters (PDEs)
+        C_M = 1.0                      # capacitance
+        temperature = 300e3            # temperature (m K)
+        R = 8.314e3                    # Gas Constant (m J/(K mol))
+        F = 96485e3                    # Faraday's constant (mC/ mol)
+
+        D_Na = Constant(1.33e-8)       # diffusion coefficients Na (cm/ms)
+        D_K = Constant(1.96e-8)        # diffusion coefficients K (cm/ms)
+        D_Cl = Constant(2.03e-8)       # diffusion coefficients Cl (cm/ms)
+        psi = F / (R * temperature)    # shorthand
+        C_phi = C_M / dt               # shorthand
+
+        Na_e_init = Constant(99.99996207620596)
+        K_e_init = Constant(4.000072443420079)
+        Cl_e_init = Constant(104.00003451962611)
+        Na_i_init = Constant(12.000067694452486)
+        K_i_init = Constant(124.9998827540323)
+        Cl_i_init = Constant(136.9999504484848)
+        phi_M_init = Constant(-68.31807655880984)
+
+        # global kappa
+        #psi = F/(R * temperature)
+        #kappa_i = F * psi * (float(D_Na) * float(Na_i_init) + float(D_K) * float(K_i_init) * float(D_Cl) * float(Cl_i_init))
+        #kappa_e = F * psi * (float(D_Na) * float(Na_e_init) + float(D_K) * float(K_e_init) * float(D_Cl) * float(Cl_e_init))
+        #E_K = float((R * temperature)/F)
+
+        #print("kappa_i", kappa_i)
+
+        #print(kappa_e)
+        #print(E_K)
+        #print(float((F*F)/(R*temperature)))
+
+        #sys.exit(0)
+
+        # Set physical parameters
+        params = namedtuple('params', ('dt', 'n_steps_ODE', 'F', 'psi', \
+                    'phi_M_init', 'C_phi', 'C_M', 'R', 'temperature'))(dt, \
+                    n_steps_ODE, F, psi, phi_M_init, C_phi, C_M, R, temperature)
+
+        # diffusion coefficients for each sub-domain
+        D_Na_sub = {1:D_Na, 0:D_Na}
+        D_K_sub = {1:D_K, 0:D_K}
+        D_Cl_sub = {1:D_Cl, 0:D_Cl}
+
+        # initial concentrations for each sub-domain
+        Na_init_sub = {1:Na_i_init, 0:Na_e_init}
+        K_init_sub = {1:K_i_init, 0:K_e_init}
+        Cl_init_sub = {1:Cl_i_init, 0:Cl_e_init}
+
+        # Create ions (channel conductivity is set below in the membrane model)
+        Na = {'c_init_sub':Na_init_sub, 'bdry': Constant((0, 0)),
+              'z':1.0, 'name':'Na', 'D_sub':D_Na_sub}
+        K = {'c_init_sub':K_init_sub, 'bdry': Constant((0, 0)),
+             'z':1.0, 'name':'K', 'D_sub':D_K_sub}
+        Cl = {'c_init_sub':Cl_init_sub, 'bdry': Constant((0, 0)),
+              'z':-1.0, 'name':'Cl', 'D_sub':D_Cl_sub}
+
+        # Create ion list. NB! The last ion in list will be eliminated, and
+        # should be the ion with the smallest diffusion coefficient D_k
+        ion_list = [K, Cl, Na]
+
+        # Membrane parameters
+        g_syn_bar = 10                       # synaptic conductivity (mS/cm**2)
+
+        # Set stimulus ODE
+        stimulus = {'stim_amplitude': g_syn_bar}
+        stimulus_locator = lambda x: (x[0] < 20.0e-4)
+
+        stim_params = namedtuple('membrane_params', ('g_syn_bar', \
+                                 'stimulus', 'stimulus_locator'))(g_syn_bar, \
+                                  stimulus, stimulus_locator)
+
+        # Get mesh, subdomains, surfaces paths
+        mesh_prefix = "meshes/3D/"
+        mesh_path = mesh_prefix + "mesh_" + str(resolution) + ".xml"
+        subdomains_path = mesh_prefix + "subdomains_" + str(resolution) + ".xml"
+        surfaces_path = mesh_prefix + "surfaces_" + str(resolution) + ".xml"
+
+        # Generate mesh if it does not exist
+        if not os.path.isfile(mesh_path):
+            from make_mesh_3D import main
+            main(["-r", str(resolution), "-d", mesh_prefix])
+
+        mesh = Mesh(mesh_path)
+        subdomains = MeshFunction('size_t', mesh, subdomains_path)
+        surfaces = MeshFunction('size_t', mesh, surfaces_path)
+
+        # Set solver parameters EMI (True is direct, and False is iterate) 
+        direct_emi = False
+        rtol_emi = 1E-5
+        atol_emi = 1E-40
+        threshold_emi = 0.9
+
+        # Set solver parameters KNP (True is direct, and False is iterate)
+        direct_knp = False
+        rtol_knp = 1E-7
+        atol_knp = 2E-40
+        threshold_knp = 0.75
+
+        # Set parameters
+        solver_params = namedtuple('solver_params', ('direct_emi',
+                                   'direct_knp', 'resolution',
+                                   'rtol_emi', 'rtol_knp',
+                                   'atol_emi', 'atol_knp',
+                                   'threshold_emi', 'threshold_knp'
+                                   ))(direct_emi, direct_knp, resolution, \
+                                      rtol_emi, rtol_knp, atol_emi, atol_knp, \
+                                      threshold_emi, threshold_knp)
+
+        # Set solver parameters (True is direct, and False is iterate)
+        #direct_emi = False
+        #direct_knp = False
+
+        # Set solver parameters
+        #solver_params = namedtuple('solver_params', ('direct_emi',
+            #'direct_knp', 'resolution'))(direct_emi, direct_knp, resolution)
+
+        # File for results
+        fname = "results/data/3D_emi/"
+
+        # Dictionary with membrane models (key is facet tag, value is ode model)
+        ode_models = {1: mm_hh}
+
+        # Solve system
+        S = SolverEMI(params, ion_list)                    # create solver
+        S.setup_domain(mesh, subdomains, surfaces)      # setup meshes
+        S.setup_parameters()                            # setup physical parameters
+        S.setup_FEM_spaces()                            # setup function spaces and numerical parameters
+        S.setup_membrane_model(stim_params, ode_models) # setup membrane model(s)
+        S.solve_system_active(Tstop, t, solver_params, filename=fname) # solve
diff --git a/elec-prop-test/unit_test.py b/elec-prop-test/unit_test.py
new file mode 100644
index 0000000..dad3672
--- /dev/null
+++ b/elec-prop-test/unit_test.py
@@ -0,0 +1,6 @@
+from astropy import units as u
+
+SI = u.S/u.m
+new = u.mS/u.cm
+
+print(SI.to(new, 1))