function [a,b,a_trace,aM] = RD_linear(ncells,converge_eps,max_steps,...
   ainit,binit,diff1,diff2,beta_rand,react_speed)
%
% [a,b,a_trace,aM]=RD_linear(60,0.0001,1200,4,4,0.25,0.0625,0.05,1);

a = ainit*ones(ncells,1);
a_trace = a';
b = binit*ones(ncells,1);
beta_init = 12;
da = zeros(ncells,1);
db = zeros(ncells,1);
beta = beta_init + (-beta_rand+2*beta_rand*rand(ncells,1));

rsp = react_speed/16;
step = 0;
err = converge_eps+1;

wb = waitbar(0,'RD Process');
while (step<max_steps) & (err>converge_eps)
    waitbar(step/max_steps);
   step = step+1;
   for c = 1:ncells
      c0 = rem(((c-1)+ncells-1),ncells)+1;
      c1 = rem(((c-1)+1),ncells)+1;
      adiff = a(c0) + a(c1) - 2*a(c);
      bdiff = b(c0) + b(c1) - 2*b(c);
      da(c) = rsp*(16-a(c)*b(c))+adiff*diff1;
      db(c) = rsp*(a(c)*b(c)-b(c)-beta(c))+bdiff*diff2;
   end
   for c = 1:ncells
      anext = a(c)+da(c);
      bnext = b(c)+db(c);
      if bnext<0
         bnext = 0;
      end
      if a(c)~=0
         erra(c) = abs(anext-a(c))/a(c);
      else
         erra(c) = anext;
      end
      if b(c)~=0
         errb(c) = abs(bnext-b(c))/b(c);
      else
         errb(c) = bnext;
      end
      a(c) = anext;
      b(c) = bnext;
   end
   a_trace = [a_trace; a'];
   clf
   plot(0,0,'k.');
   hold on
   xlabel('Cell Locations');
   ylabel('Morphogen Concentration');
   plot(ncells+1,10,'k.');
   plot(a);
   aM(step) = getframe;
   err = max(max(erra),max(errb));
end
close(wb);
