solving a second spatial derivative using ode45

I am trying to solve this above system using ode45 but I am confused on how to code for solving the second spatial derivative. The way this problem is set up is a 1D cabel model with a certain set of nodes, spacing between the nodes, and specific nodes where a stimulus is applied. So the above ODE's are being solved at each node.

function dUdt = fhn(t,U,stim,numNodes,dx,stimNodes,u0)

% FitzHugh-Nagumo equations defined to be plugged into ode45

%

% INPUT: t time

% U vector that holds u and v

% numNodes number of nodes along the fiber

% dx spatial step, spacing between nodes

% stimNodes set nodes where stim is applied

%

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

% v at each node

% Other needed variables

a = 0.13;

b = 0.013;

c1 = 0.26;

c2 = 0.1;

D = 5;

% There are 2 ODE's per node

% Using no-flux boundary conditions on the two ends of the cable

% At the first and last nodes, use a special version of the diffusion term

% in which du/dx is 0

for i = 1:numNodes

% If at the first node solve the system using no flux boundary conditions

if i == 1

dudt = (c1*U(1))*(U(1)-a)*(1-U(1))-(c2*U(1)*U(2))+D*((U(1)-u0)/dx^2);

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

% If at the last node then solve with no flux boundary conditions

elseif i == numNodes

dudt = (c1*U(1))*(U(1)-a)*(1-U(1))-(c2*U(1)*U(2))+D*((U(1)-u0)/dx^2);

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

% If at one of the specified stimulus nodes then the stim term needs to

% be incorporated

elseif i == stimNodes

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

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

% Otherwise just solve the system normally without stimulus or no flux

% boundary conditions

else

dudt = (c1*U(1))*(U(1)-a)*(1-U(1))-(c2*U(1)*U(2))+(U(1(i+1))/dx^2);

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

end

end

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

% v

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

dUdt = [dudt; dvdt];

Above is my code so far, but I am confused on how to solve the second spatial derivative since the equation given to solve for it on a 1D cable model uses central finite differences. The equation is ( u(n+1) + u(n-1) - 2u ) / dx^2 where u corresponds to U(1) in the code and n corresponds to the number of nodes.

Answer 1

[J. Alex Lee]

This comment is incorrect:

% Create a two element vector that holds the derivative equations of u and
% v
% U(1) is u and U(2) is v
dUdt = [dudt; dvdt];

You say: "So the above ODE's are being solved at each node", but that's not quite right either. You are discretizing your spatial derivatives into finite differences (you're wanting to implement a method of lines), so that you have numNodes equations for the discretized , but is still a scalar. So in total, your fhn() function should return a vector of length numNodes+1.

Answer 2

[darova]

You difference scheme is wrong

dudt = (c1*U(1))*(U(1)-a)*(1-U(1))-(c2*U(1)*U(2))+D*((U(1)-u0)/dx^2)

Here is correct version

uu = u(i,j);
vv = v(i,j)
(u(i+1,j)-uu)/dt = c1*uu*(uu-a)*(1-uu) - c2*uu*vv + stim + D*(u(i,j+1)-2*uu+u(i,j-1))/dx^2
u(i+1,j) = LONG_EXPRESSION*dt + uu;

Then for loop will

for i = 1:m-1
    for j = 2:n-1
        u(i+1,j) = LONG_EXPRESSION*dt + u(i,j);
        v(i+1,j) = b*(u(i,j)-v(i,j))*dt + v(i,j); 
    end
end

u and v - matrices

Do you have more initial/boundary conditions except of this?

% Using no-flux boundary conditions on the two ends of the cable
% At the first and last nodes, use a special version of the diffusion term
% in which du/dx is 0