function [ v ] =nma_solve_reaction_ODE( a,v_old,w,delt, tolerance)
%solves the reaction ODE part of the FitzHugh-Nagumo equations
%nma_solve_reaction_ODE solves the reaction ODE part of the
%FitzHugh-Nagumo equations.
%
%Backward Euler is used in conjugtion of Newton root finding to solve
%for the solution of the implicit equation generated from the backward
%Euler scheme.
%
%INPUT
% a: problem parameter, a scalar
% v_old: initial data, 2D grid data
% w: 2D data, represent the w variable
% delt: the time step to use
% tolerance: a number, used by Newton root finder to know when to stop.
%
%OUTPUT
% v: The value of v after integrating the ODE for delt time
%
%HW2, problem 3, Math 228B, UC Davis
%by Nasser M. Abbasi
n = size(v_old,1);
I = 0; %current is zero in this problem.
M = v_old-delt*w+delt*I;
v = M; %initial guess for Newton root finder
MAX = 10; %upper limit on how many time to look for root per grid point
%should never be reached. Message will be displayed if so.
%LOOP over each grid point, and integrate the equation for that point
for i = 1:n
for j = 1:n
k = 0; %count how many times Newton loop run for
done = false;
while not(done)
num = v(i,j) - delt*(a-v(i,j))*(v(i,j)-1)*v(i,j);
den = 1+delt*v(i,j)*(3*v(i,j)-2)*v(i,j)+a*(delt-2*delt*v(i,j));
v_new = v(i,j) - (num-M(i,j))/den;
if abs(v_new-v(i,j)) < tolerance
done = true;
end
v(i,j) = v_new; %update root for next Newton iteration
%check if we looped too many times. If so, something wrong in code
k = k+1;
if k>MAX
done = true;
fprintf('WARNING, root not found\n');
end
end
end
end
end