%% FitzHugh-Nagumo avec diffusion 1D

% parametre des equations
epsilon = 0.08; % 0.08;
b = 0.7;
c = 0.8;
D = .01;
I = 0.0;

% parametres de simulation, espace
% x in [0,1]
h = 0.1;
x0 = 0;
x1 = 10;
x = x0:h:x1; 
J = length(x);

% variable dynamiques
v = zeros(J,1); % stocke seulement l'etat au temps t
w = zeros(J,1);
newv = zeros(J,1);
neww = zeros(J,1);

% discretisation du Laplacien 
neumann = 0;
periodique = 1;
dirichlet = 0;

if neumann % Condition de Neumann v' = 0
    L = sparse(1:J,1:J,-2); % matrice creuse, compacte en memoire
    L = spdiags(ones(J,2),[-1 1],L); % forme la matrice tridiagonale
    L(1,:) = 0;
    L(J,:) = 0;
elseif periodique % Condition periodiques
    L = sparse(1:J,1:J,-2); % matrice creuse, compacte en memoire
    L = spdiags(ones(J,2),[-1 1],L); % forme la matrice tridiagonale
    L(1,J) = 1;
    L(J,1) = 1;
elseif dirichlet
    
end

C = sparse(1:J,1:J,1); 
if neumann % Condition de Neumann v' = 0
    C(1,1) = 0;
    C(1,2) = 1;
    C(J,J) = 0;
    C(J,J-1) = 1;
elseif periodique % Conditions periodiques
% rien a faire
end

% condition initiales
if neumann 
    v(x <= (0.2*x1) ) = 0.5;
    v(x > (0.2*x1) )  = -1.1;
    w(1:J) = -0.5;    
elseif periodique
    v(x <= (0.2*x1) & x >= (0.15*x1) ) = 1.0;
    v(x > (0.2*x1) | x < (0.15*x1) )  = -1.1;
    w(x <= (0.15*x1)) = 0.0;
    w(x > (0.15*x1)) = -0.5;
elseif dirichlet
    % --
end

% parametres de simulation, temps
t0 = 0;
tfinal = 200; 
t = t0;
% k doit etre < a h^2/2/D
k = min(0.5,0.5*( h^2/2/D ));
 

figure(1); clf;
plot(x,v,x,w);
axis([x0 x1 -2 2])
pause

% BOUCLE PRINCIPALE
while t < tfinal
    drawnow;
    newv = C*v + k*(v-v.^3/3-w+I) + ...
        D*k/h^2*L*v;
    neww = w + k*epsilon*(v + b - c*w);    
    v = newv;
    w = neww;
    plot(x,v,x,w);
    axis([x0 x1 -2 2])
    t = t + k
end