function [a,b] = Turingk(steps,ncells)
%

react_speed = 1.0;
rsp = react_speed/16;
ainit = 4;
binit = 4;
beta_init = 12;
beta_rand = 0.05;
diff1 = 0.25;
diff2 = 0.0625;
a = zeros(1,ncells);
b = zeros(1,ncells);
da = a;
db = b;
beta = a;

for index = 1:ncells
   a(index) = ainit;
   b(index) = binit;
   beta(index) = beta_init + frand(-beta_rand,beta_rand);
end

for s = 1:steps
   for nc = 0:ncells-1
      i0 = rem((nc+ncells-1),ncells)+1;
      i1 = rem((nc+1),ncells)+1;
      nm = nc+1;
      adiff = a(i0) + a(i1) - 2*a(nm);
      bdiff = b(i0) + b(i1) - 2*b(nm);
      da(nm) = rsp*(16-a(nm)*b(nm)) + adiff*diff1;
      db(nm) = rsp*(a(nm)*b(nm) - b(nm) - beta(nm)) + bdiff*diff2;
   end
   for nc = 1:ncells
      a(nc) = a(nc) + da(nc);
      b(nc) = b(nc) + db(nc);
      if b(nc)<0
         b(nc) = 0;
      end
   end
end
