Nonlinear 2D modeling of Neural Electromagnetic field in Matlab

I am trying to replicate the MATLAB simulation presented in this paper. More specifically, I have to code the solution to this equation

$$\frac{a}{2R_i}\frac{\partial^2 V_m}{\partial x^2} - \left[g_\text{Na}m^3 h(V_m - E_\text{Na}) + g_\text{K} n^4(V_m - E_\text{K}) + g_\text{S}(V_m - E_\text{S})\right] = C_m \frac{\partial V_m}{\partial t} + \frac{a}{2R_i}\frac{\partial E_x}{\partial x}(x,t) \enspace .$$

Here $m$, $n$, and $h$ are gating variables that are ordinary differential equations and are functions of time, and I have solved for them already. I have also written the difference equation for this PDE. I would like to set $\frac{\partial E}{\partial x}(x,t)$ equal to 1 (denoted as $g$ in the script) for now (and tackle that next). I am using a finite differences method (I think)

clear all;
T=1; % Max time (s)
L=1; % Max length of neuron
dt=0.04; %time step
dx=0.01; %space step
t=[0:dt:T]; %time vector
x=[-L:dx:L]; %Space vector

M=length(x);
g=1; %The dE/dx bit
R_i=1; %Intracellular resistance
A=a./2*R_i;
I=0.1; %External Current Applied
ENa=55.17; % mv Na reversal potential
EK=-72.14; % mv K reversal potential
El=-49.42; % mv Leakage reversal potential
g_Na=1.2; % mS/cm^2 Na conductance
g_K=0.36; % mS/cm^2 K conductance
g_l=0.003; % mS/cm^2 Leakage conductance
b=4;
Cm=0.01; %Membrane capacitance
C=Cm;
%Solve for m,n,h

%Solve for m, n, h
V=-60; % Initial Membrane voltage
m=am(V)/(am(V)+bm(V)); % Initial m-value
n=an(V)/(an(V)+bn(V)); % Initial n-value
h=ah(V)/(ah(V)+bh(V)); % Initial h-value
y0=[V;n;m;h];
tspan = [0,max(t)];
%Matlab's ode45 function
[time,V] = ode45(@HH,tspan,y0);
V_hh=V(:,1);
n=V(:,2);
m=V(:,3);
h=V(:,4);
N=length(time);
plot(time,m,time,h,time,n);
V_initialgrid=zeros(M,N);
figure();
plot(time,V_hh)
V_initialgrid(1,:)=V_hh';

xlabel('time');
ylabel('length');
for j=2:N-1 %Stepping through time
B=g_Na*m(j).^3*h(j)+g_K*n(j).^4+g_l;
for i=2; %stepping thru space
V_initialgrid(i,j)=(A* V_initialgrid(i-1,j)+B* V_initialgrid(i,j)+C* V_initialgrid(i-1,j)+g);

for i=3:M-1 %stepping thru space
V_initialgrid(i+1,j)=(-dt*(g-(A.*(( V_initialgrid(i,j-1)-2.* V_initialgrid(i,j)+ V_initialgrid(i,j+1))./(dx.^2))-B.*-60))./C)+ V_initialgrid(i,j);
end
end
end
imagesc(log(abs(V_initialgrid)))

xlabel('Time');
ylabel('L');


The boundary conditions are as follows:

$$\frac{\partial V}{\partial x}=\frac{\partial m}{\partial x}=\frac{\partial h}{\partial x}=\frac{\partial n}{\partial x}=0\quad \text{for } x=\pm L$$

and

$$\frac{\partial V}{\partial t}=\frac{\partial m}{\partial t}=\frac{\partial h}{\partial t}=\frac{\partial n}{\partial t}=0 \enspace .$$

I know my error lies within the indexing of the solution of the difference equation.

• Welcome to SciComp.SE. Try to include the information that is necessary to address your question. I just included the equation you want to solve, please check that is correct. – nicoguaro Jun 22 '15 at 19:42
• Since you consider $m,n,h,E$ to be given, this is really just the heat equation. There is a nice and detailed explanation of how to solve it, e.g. in Chapter 2 of LeVeque's book on finite differences. – David Ketcheson Jun 23 '15 at 8:13