solving 1D cable model with ode45

I'm confused on how to update my 0-dimensional fitzhugh-nagumo code to solve the system above with a 1D cable model. For the zero-dimensional case, you only had to solve 2 ODEs. For this problem, there are two ODEs per node. You're supposed to use no-flux boundary conditions on the two ends of the cable. In other words, at the first and last nodes, use a special version of the diffusion term in which du/dx is 0. Input should be updated to specify the number and spacing of nodes as well as the set nodes at which the stimulus is applied. The program should output a time vector and separate arrays for u and v. The u and v arrays should have one row per timestep, and one column per node on the fiber.

Below is my 0D code.

function [t,u,v] = fhn0D(tf,u0,v0,t0,stim,tstim)

% Solve the FitzHugh-Nagumo equations for a single cell (i.e., without

% spatial coupling)

% du/dt = (c1*u)*(u - a)*(1 - u)*(c2*u*v) + stim

% dv/dt = b(u - v)

% stim is a stimulus current that can be applied for a short time at

% the beginning of the simulation

% u represents membrane potential and ranges from 0 (rest) to 1 (excited)

% v is a recovery variable in the same range (also in 0-1 range)

% t is time in milliseconds

% INPUT: tf duration of simulation

% u0 initial value of u

% v0 initial value of v

% t0 initial value of t

% stim strength of stimulus

% tstim time of applied stimulus

%

% OUTPUT: t time vector (ms)

% u membrane potential vector

% v recovery vector

[t,U] = ode45(@(t,U)fhn(t,U,stim,tstim),[t0 tf],[u0; v0]);

% First and second column of U correspond to u and v respectively

u = U(:,1);

v = U(:,2);

% Plot u and v vs. t

plot(t,u,t,v)

title('Solution of FitzHugh-Nagumo Equations for a Single Cell with ODE45');

xlabel('Time t');

ylabel('Solution U');

legend('u','v');

end

function dUdt = fhn(t,U,stim,tstim)

% FitzHugh-Nagumo equations defined to be plugged into ode45

%

% INPUT: t time

% U vector that holds u and v

%

% OUTPUT: dUdt vector containing solutions to derivatives of u and v

% Other needed variables

a = 0.13;

b = 0.013;

c1 = 0.26;

c2 = 0.1;

% Conductivity constant

D = 5;

% Preallocate the output vector

dUdt = zeros(2,1);

% Only want the stimulus applied over a certain time point of the action

% potential not for the entire action potential

% If time is in the range of tstim then stim is incorporated into the

% equation otherwise stim is 0

if t <= tstim

dudt = (c1*U(1))*(U(1)-a)*(1-U(1))-(c2*U(1)*U(2))+ stim;

else

dudt = (c1*U(1))*(U(1)-a)*(1-U(1))-(c2*U(1)*U(2));

end

% Create a two element vector that holds the derivative equations of u and

% v

% U(1) is u and U(2) is v

dvdt = b*(U(1)-U(2));

dUdt = [dudt; dvdt];

end