D_a = 0.01;
D_b = 0.3;
D_c = 0.002;

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

b_a = 0.01;
b_b = 0.03;
b_c = 0;

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

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

KT = 450;
KP = 6;
KY = 80;

dt = 1;
h = 1;


% FIXME: this works only because G_* is 1
a = ones(KT, KY);
b = ones(KT, KY);
c = ones(KT, KY);



dra = 1 - r_a - 2 * D_a;
drb = 1 - r_b - 2 * D_b;
drc = 1 - r_c - 2 * D_c;

na = zeros(1, KY);
nb = zeros(1, KY);
nc = zeros(1, KY);
        
s = r_a * (1 + 0.1 * (rand(1,KY)-0.5));

for n=1:KT
    for k=1:KP
        oa = a(n,:);
        ob = b(n,:);
        oc = c(n,:);
        
        for j=1:KY
            ooa = oa(j);
            oob = ob(j);
            ooc = oc(j);
            
            if j == 1
                j_links = j;
            else
                j_links = j - 1;
            end
            if j == KY
                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 = dra * ooa + D_a * (oa(j_links) + oa(j_rechts));
            olddecaydiffB = drb * oob + D_b * (ob(j_links) + ob(j_rechts));
            olddecaydiffC = drc * ooc + D_c * (oc(j_links) + oc(j_rechts));
            na(j) = olddecaydiffA + aq;
            nb(j) = olddecaydiffB - aq + b_b;
            nc(j) = olddecaydiffC + r_c * ooa;
        end
        a(n,:) = na;
        b(n,:) = nb;
        c(n,:) = nc; 
    end
    a(n+1,:) = na;
    b(n+1,:) = nb;
    c(n+1,:) = nc; 
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)])