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

T = KT*KP;
dt = T / steps;
L = 1;
N = 80;
h = 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;


% FIXME: this works only because G_* is 1
a = ones(steps/KP, N);
b = ones(steps/KP, N);
c = ones(steps/KP, N);

na = zeros(1, N);
nb = zeros(1, N);
nc = zeros(1, N);
        
%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);

for n=1:steps/KP
    for k=1:KP
        oa = a(n,:);
        ob = b(n,:);
        oc = c(n,:);
        
        if 0        
        for j=1:N
            ooa = oa(j);
            oob = ob(j);
            ooc = oc(j);
            
            if j == 1
                j_links = j;
            else
                j_links = j - 1;
            end
            if j == N
                j_rechts = j;
            else
                j_rechts = j + 1;
            end
            
            aq = s(j) * oob * ((ooa^2 + b_a) / (1 + s_a * ooa^2)) / (s_b + s_c * ooc);
            olddecaydiffA = -r_a * ooa + D_a / h^2 * (oa(j_links) + oa(j_rechts) - 2 * ooa);
            olddecaydiffB = -r_b * oob + D_b / h^2 * (ob(j_links) + ob(j_rechts) - 2 * oob);
            olddecaydiffC = -r_c * ooc + D_c / h^2 * (oc(j_links) + oc(j_rechts) - 2 * ooc);
            na(j) = olddecaydiffA + aq;
            nb(j) = olddecaydiffB - aq + b_b;
            nc(j) = olddecaydiffC + r_c * ooa;
        end
        else
        aq = s .* ob .* ((oa.^2 + b_a) ./ (1 + s_a * oa.^2)) ./ (s_b + s_c * oc);
        olddecaydiffA = -r_a * oa + D_a * oa * Lap;
        olddecaydiffB = -r_b * ob + D_b * ob * Lap;
        olddecaydiffC = -r_c * oc + D_c * oc * Lap;
        
        na = olddecaydiffA + aq;
        nb = olddecaydiffB - aq + b_b;
        nc = olddecaydiffC + r_c * oa;
        end
        
        a(n,:) = oa + dt * na;
        b(n,:) = ob + dt * nb;
        c(n,:) = oc + dt * nc;
    end
    a(n+1,:) = a(n,:);
    b(n+1,:) = b(n,:);
    c(n+1,:) = c(n,:); 
end

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)])