% 1D time dependent activator-inhibitor system on periodic domain

L = 4;		% length of periodic domain
N = 256;	% space sample points
LifeT = 1;	% life time of the sea shell
nT = 256;	% time sample points

% system parameters
r_a = 0.02;	% activator removal rate
D_a = 0.01;	% activator diffusion coefficient
b_a = 0.001;	% basic activator production
r_b = 0.03;	% inhibitor removal rate
D_b = 0.4;	% inhibitor diffusion coefficient
b_b = 0.0;	% basic inhibitor production
s = r_a * (0.99 + 0.2 * rand(1,N));	% activator production rate

% help variables
h = L/N;	% grid size
dt = LifeT/nT;	% time step


% Build Laplace Operator (Finite Differences)
e = ones(N,1) / h^2;
Lap = spdiags([e -2*e e], -1:1, N, N);
Lap(N,1) = e(1);
Lap(1,N) = e(1);

% set initial values
a = zeros(nT, N);
b = zeros(nT, N);
a(1,:) = 1.0; % + 0.01 * rand(1,N);
b(1,:) = 1.0; % + 0.01 * rand(1,N);

% solve problem (Fixed Step Euler Forward)
for n=1:nT-1
  oa = a(n,:);
  ob = b(n,:);
  a(n+1,:) = oa + dt * (s.*(oa.^2./ob + b_a) - r_a.*oa + D_a*oa*Lap);
  b(n+1,:) = ob + dt * (s.*oa.^2 - r_b*ob + D_b*ob*Lap + b_b);
end

colormap([[0:0.5/63:0.5]',[0:0.7/63:0.7]',[0:1/63:1]'])
subplot(1,2,1), image(a);
subplot(1,2,2), image(b);