global N s D_a D_b D_c b_a b_b b_c s_a s_b s_c r_a r_b r_c Lap indA indB indC

KT = 200;
KP = 6;
KY = 80;
steps = KT*KP;

T = KT*KP;
dt = T / steps;
L = 1;
N = 400;
h = 1/N;

indA = 1:N;
indB = N + (1:N);
indC = 2*N + (1:N);

D_a = 0.01/KY^2;
D_b = 0.3/KY^2;
D_c = 0.002/KY^2;

r_a = 0.08;
r_b = 0;
r_c = 0.01;

b_a = 0.01;
b_b = 0.06;
b_c = 0;

s_a = 0.1;
s_b = 0;
s_c = 1;

G_a = 1;
G_b = 1;
G_c = 1;


s = r_a * (1 + 0.1 * (rand(1,N)-0.5));
%s = r_a * (1 + 0.05 * sin(linspace(0, 10*pi, N)))';

s = s';

e = ones(N,1)/h^2;
Lap = spdiags([e -2*e e], -1:1, N, N);
Lap(1,1) = -e(1);
Lap(N,N) = -e(1);

y0 = ones(3*N,1);

opt = odeset('JPattern', Lap, 'RelTol', 1e-2);
sol = ode15s(@matlab_muschel_f, [0, T], y0);

t = linspace(0, T, 200);
y = deval(sol, t);
a = y(indA,:);
b = y(indB,:);
c = y(indC,:);

%subplot(1,3,1), imagesc(a);
%subplot(1,3,2), imagesc(b);
%subplot(1,3,3), imagesc(c);

imagesc(a');
color = [makecolor(135/255), makecolor(60/255), makecolor(20/255)];

colormap(color)