% test it does correct thing Schnakenberg equation
%   https://thatsmaths.com/2013/04/25/spots-and-stripes/
%   see also schnakenberg2d.m for parameters, from http://homepage.univie.ac.at/marie-therese.wolfram/teaching.html
%   slightly adapted to match parameters for 1D
%   
clear
seed = 3
rand('state', seed)

num_x = 60;      % number of grid points in x
num_t = 10000;   % number of time steps
dx = 1/num_x;
x = 0: dx: 1-dx;
dt = 0.00001;   

D = diag([1.0, 4.8]);
gamma = 1000; % or 100, 1600, 6400
a = -0.55;
b  = 1.9;
reaction_f = @(r) gamma*[ a - r(1) + r(1)^2*r(2); 
                          b        - r(1)^2*r(2)]; 
sigma = 0.01;
u0(1,:) = (a+b) +     sigma*rand(1,num_x);
u0(2,:) = b/(a+b)^2 + sigma*rand(1,num_x); 
u = reaction_diffusion( reaction_f, D, u0, num_t, dt);
  
figure(1)
clf()
hold on
p1 = plot(x, u0, '--', 'linewidth', 3);
p2 = plot(x, u, 'linewidth', 3);
set(gca, 'fontsize', 18, 'linewidth', 3)
ell = legend([p1; p2], 'initial $u_1$', 'initial $u_2$', 'final $u_1$', 'final $u_2$', 'Location', 'EastOutside')
set(ell,'Interpreter','Latex')
ylabel('reagent concentrations');
xlabel('x');

set(gcf,'PaperUnits', 'centimeters', ...
	'PaperPosition', [0 0 20 12], ...
	'PaperSize', [20 12]);
filename = sprintf('Plots/reaction_diffusion_test_4_%03d.pdf', seed);
print('-dpdf', filename);
