function nma_HW2_math_228B_problem3()
% solves the FitzHugh-Nagumo on unit square
% Math 228B, UC Davis, winter 2011.
%
% This function solves the FitzHugh-Nagumo equations on unit
% square
% V_t = D ( V_xx+V_yy) + (a-V)(v-1)v - w + I
% w_t = epsilon(V- gamma*w)
%
% subject to Nuemann homogenous boundary conditions.
% The solver splitting method to solve V_t = D ( V_xx+V_yy)
% and then solve V_t = (a-V)(v-1)v - w + I. The ADI scheme is
% used to solve V_t = D ( V_xx+V_yy), backward Euler is used to
% solved the reaction equation V_t = (a-V)(v-1)v - w + I.
% To solve w_t = epsilon(V- gamma*w), forward Euler is used.
%
% See my HW2 report for more information.
% by Nasser M. Abbasi 2/14/2011
%PARAMETERS
a = 0.1;
gamma = 2;
epsilon = 0.005;
D = 5*10^-5; %Diffusion constant
h = 0.01; %space step
I = 0; %zero current
MODE = 1; %set this to 2 to plot the W state variable. 1 to plot V .
GENERATE_GIF=false; %turn to true to generate animated gif
% Select the part(b) or part(c) initial data by commenting
% the part needed
% part(b)
ic_v = @(X,Y) exp(-100*(X.^2+Y.^2));
ic_w = @(X,Y) zeros(size(X));
%part(c)
%ic_v = @(X,Y) 1-2*X;
%ic_w = @(X,Y) 0.05*Y;
% Allocate mesh for plotting
[X,Y] = meshgrid(h/2:h:1-h/2,h/2:h:1-h/2);
% generate initial data for w and voltage V
w = ic_w(X,Y);
v = ic_v(X,Y);
% These below are for plotting only to find good limit for z axis
maxv = 1.5*max(max(v));
minv = -0.5+min(min(v));
% Select time step. Use same as space step to keep second order
% since splitting method is used.
delt = h;
% Select a tolerance to use for Newton root solver to use for solving
% the reaction ODE
tolerance = 10^-7;
% Find number of steps to run the simulation based on number of seconds
number_of_seconds = 300; % use 600 for part (d)
number_of_steps = round(number_of_seconds/delt);
% Generate the A,B Matrices to use for ADI solver. THis is done once
% and not changed, since the same grid size is used for each step.
[A,A_rhs]= ...
nma_generate_A_and_ARHS_for_2D_diffusion_Neumman(round(1/h),D,delt,h);
% Get ready to start main loop. Initialize counter and set the time step
stepA = delt;
step_number = 0;
fig_handle = figure();
set(fig_handle,'Name',....
'FitzHugh-Nagumo equations simulation. HW2,Math 228B, Nasser M. Abbasi');
frame_number = 0;
[im,map] = initialize_plot(fig_handle,step_number,X,Y,v,minv,maxv,delt,GENERATE_GIF);
while step_number