clear all; close all; clc;

% interactions:
M  = [-2  -3;
       1   1];
% diffusion:
D  = [8 1]*0.1;
% number of cells:
N  = 200;
% distance between cells:
dx = 0.1;

% det M (if positive, uniform solution could be stable)
det(M)
% det M at max  (if negative, there is a wavelength at which system becomes unstable)
det(M) - (M(1,1)*D(2) + M(2,2)*D(1))^2/D(1)/D(2)/4
% wavelength
sqrt(2 * D(1) * D(2) / (D(1) * M(2,2) + D(2) * M(1,1))) * 2 * pi

% plot larger eigenvalue
pluseig=[];
for k2 = 0:0.1:4
    pluseig = [pluseig max(real(eig(M - k2*diag(D))))];
end

% plot determinant
dets=[];
for k2 = 0:0.1:4
    dets = [dets det(M - k2*diag(D))];
end

% time span to simulate system over
tspan = [0 500];

% start with random positive numbers for the concentrations
y0 = 0.01*randn(N*2,1);
y0 = y0 - min(y0);


[t,y] = ode45(@(t,y) Turing_ode(t,y,M,D,N,dx), tspan, y0);


maxy = max(y(end,:));
miny = min(y(end,:));


figure('Position',[100 100 900 900])
hold on;
plot(t,y(:,1:N)*0+50*(1:N),'k:','LineWidth',2)
plot(t,y(:,1:N)+50*(1:N),'LineWidth',5)
% plot(t,y(1,1:4),'LineWidth',5)
hold off;


figure('Position',[100 1000 900 900])
hold on;
plot(pluseig);
hold off;


figure('Position',[1000 1000 900 900])
hold on;
plot(dets);
hold off;


figure('Position',[1000 100 900 900])
hold on;
plot((1:N)*dx,y(end,1:N))
hold off;