clc
close all hidden
clear all

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

dtau = 0.001;
tau_max = 20;

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

% beta = sech(z);
beta = exp(-10*z.^2);
% beta = exp(-z.^2/100).*(2 + cos(z))/3;
alpha = ones(size(beta));

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

hf=figure('Position',[100,100,600,600]);
rect = get(hf,'Position');
rect(1:2) = [0 0]; % these commands ensure we capture the actual figure.


plot(z,alpha,z,beta,'LineWidth',2)
    title(['t = ',num2str(0,'%2.1f')])
    axis([-50,50,0,1.2])
    ylabel('\alpha, \beta')
    xlabel('x')
    pause(0.01)
    gif_add_frame(gcf,'reaction_diffusion.gif',1);
    
% main loop, over each timestep
for n = 1:tau_max/dtau
   
    alpha = ifft(exp(-k.^2*dtau).*fft(alpha));
    beta = ifft(exp(-k.^2*dtau).*fft(beta));

    alpha = alpha.*exp(-beta*dtau);
    beta = beta.*exp(alpha*dtau);
    % into the real difft(exp(-k.^2*dz).*fft(u));omain
    if mod(n,500) == 0
    plot(z,alpha,z,beta,'LineWidth',2)
    title(['t = ',num2str(n*dtau,'%2.1f')])
    axis([-50,50,0,1.2])
    ylabel('\alpha, \beta')
    xlabel('x')
    pause(0.01)
    gif_add_frame(gcf,'reaction_diffusion.gif',1);
    end
end