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