# Gray-Scott reaction-diffusion model

clear *
gs_F = 0.026;   # Gray-Scott parameters
gs_k = 0.055;
gs_Du = 2e-5;
gs_Dv = 7e-6;
Lx = 2;         # x-length
dx = 0.01;      # x-cell size
dt = 0.5;       # integration time step
o_intv = 50;    # output (plot) interval
max_t = 3000;   # maximum time

n_intv = o_intv / dt;
n_outp = max_t / o_intv;

function gs_plot(x, uv, t)
  plot (x, uv);
  axis([0 2 -0.05 1.15]);
  title (["t = ", sprintf("%g", t)]);
  xlabel("x");
  legend("u", "v");
  refresh;
endfunction

# initial condition
x = 0 : dx : Lx - dx;
u = sin(pi * (x - 1));
v = cos(pi * (x - 1));
gs_plot (x, [u; v], 0.);

# time integration
for j = 1 : n_outp
  for jj = 1 : n_intv
    d2udx2 = (shift(u, -1) - 2 * u + shift(u, 1)) / (dx * dx);
    d2vdx2 = (shift(v, -1) - 2 * v + shift(v, 1)) / (dx * dx);
    u_new = u + (gs_Du * d2udx2 - u .* v .* v + gs_F * (1 - u)) * dt;
    v_new = v + (gs_Dv * d2vdx2 + u .* v .* v - (gs_F + gs_k) * v) * dt;
    u = u_new;
    v = v_new;
  endfor
  gs_plot (x, [u; v], j * o_intv);
endfor