global N s D_a D_b D_c D_d D_e D_f
global b_a b_b b_c b_d b_e b_f
global c_a c_b c_c c_d c_e c_f
global r_a r_b r_c r_d r_e r_f
global s_a s_b s_c s_d s_e s_f
global Lap indA indB indC indD indE indF

KT = 396;
KP = 4;
KY = 205;
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);
indD = 3*N + (1:N);
indE = 4*N + (1:N);
indF = 5*N + (1:N);

D_a = 0.0040/KY^2;
D_b = 0.0000/KY^2;
D_c = 0.0100/KY^2;
D_d = 0.4000/KY^2;
D_e = 0.0200/KY^2;
D_f = 0.2000/KY^2;

r_a = 0.01;
r_b = 0;
r_c = 0.015;
r_d = 0.02;
r_e = 0.1;
r_f = 0.001;

b_a = 0;
b_b = 0.12;
b_c = 0.005;
b_d = 0.01;
b_e = 0.005;
b_f = 0.0015;

s_a = 0.1;
s_b = 0.3;
s_c = 0;
s_d = 0;
s_e = 0.1;
s_f = 0.0005;

c_a = 0;
c_b = 0.008;
c_c = 0;
c_d = 0;
c_e = 0.15;
c_f = 0.04; % TODO: c_f(x)

G_a = 1;
G_b = 6;
G_c = 1;
G_d = 1;
G_e = 0.0001;
G_f = 0;

%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 = [ G_a * ones(N, 1);
       G_b * ones(N, 1);
       G_c * ones(N, 1);
       G_d * ones(N, 1);
       G_e * ones(N, 1);
       G_f * ones(N, 1); ];

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

t = linspace(0, T, 600);
y = deval(sol, t);

a = y(indA,:)';
b = y(indB,:)';
c = y(indC,:)';
d = y(indD,:)';
e = y(indE,:)';
f = y(indF,:)';

subplot(3,2,1), imagesc(a);
subplot(3,2,2), imagesc(b);
subplot(3,2,3), imagesc(c);
subplot(3,2,4), imagesc(d);
subplot(3,2,5), imagesc(e);
subplot(3,2,6), imagesc(f);

%colormap([makecolor(0.1), makecolor(0.7), makecolor(1)])
colormap('default');