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 = 80;
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.08;
b_c = 0;

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

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


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

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);
f0 = muschel_f(y0, 0);

t = linspace(0, T, 100);
y = dassl(@muschel_f, y0, f0, 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);

function y = makecolor(x)
% makecolor: return a vector containing 64 steps of the actual color percentage
    if (x == 0)
        y = zeros(64,1);
    else
        y = [0:x/63:x]';
    end
end

colormap([makecolor(0.1), makecolor(0.7), makecolor(1)])