Diffusion_Reaction.m

%Modified from http://www.wikiwaves.org/Reaction-Diffusion_Systems
%Reaction Diffusion
% A + E -> 2A
% Differential Equations
% A_t = A_xx + A*E
% E_t = E_xx - E*A
% Equivalente to 
% A_t = A_xx + A*(1-A)

clc
close all hidden
clear all

N = 4*256;
L = 20;
set(gca,'FontSize',18)

dt = 0.001;
t_max = 20*2;

x = L*(2*pi/N)*(-N/2:N/2-1)';

% A = sech(x);
A = exp(-10*x.^2);
% A = exp(-x.^2/100).*(2 + cos(x))/3;
E = ones(size(A));

k = 1/L*[0:N/2-1 -N/2 -N/2+1:-1]';

filename='reaction_diffusion.gif';
delete filename

plot(x,E,x,A,'LineWidth',2)
    title(['t = ',num2str(0,'%2.1f')])
    axis([-50,50,0,1.2])
    ylabel('E, A')
    xlabel('x')
    pause(0.01)
    frame = getframe(1);
    im = frame2im(frame);
    [imind,cm] = rgb2ind(im,256);
    imwrite(imind,cm,filename,'gif', 'Loopcount',inf)
    
% main loop, over each timestep
for n = 1:t_max/dt
   
    E = ifft(exp(-k.^2*dt).*fft(E));
    A = ifft(exp(-k.^2*dt).*fft(A));

    E = E.*exp(-A*dt);
    A = A.*exp(E*dt);
    % into the real difft(exp(-k.^2*dx).*fft(u));omain
    if mod(n,500) == 0
        plot(x,E,x,A,'LineWidth',2)
        title(['t = ',num2str(n*dt,'%2.1f')])
        axis([-50,50,0,1.2])
        ylabel('E, A')
        xlabel('x')
        pause(0.01)
        frame = getframe(1);
        im = frame2im(frame);
        [imind,cm] = rgb2ind(im,256);
        imwrite(imind,cm,filename,'gif','WriteMode','append')
    end
end