% m55_grayscott.m - Gray-Scott eqs. movie
%        Explicit Euler in t, finite diffs. in x
%        Qualitatively but not quantitatively correct.

  hold off, shg

% Grid and initial condition:
  J = 60;
  h = 1/J;                                  % space step
  k = 2000*h^2;                             % time step
  s = h*(1:J)';
  tplot = 100;
  nplot = round(tplot/k);
  k = tplot/nplot;                          
  [xx,yy] = meshgrid(s,s);  
  x = xx(:); y = yy(:); 
  du = sqrt((x-.2).^2+(y-.2).^2);
  dv = sqrt((x-.3).^2+2*(y-.3).^2);         % initial
  u = 10*du; u(u>1) = 1;                    % conditions
  v = 1-10*dv; v(v<0) = 0;
  epsu = 5e-5; epsv = 2e-5; c = .065;       % constants
  F = input('F? (e.g. .03, .06, .044) ');

% Construct time-stepping matrices:
  I = speye(J); II = speye(J^2);         
  D = h^(-2)*toeplitz([-2 1 zeros(1,J-2)]); % 2nd-order derivative
  D(1,end) = D(2,1); D(end,1) = D(1,2);     % periodic BCs
  L = kron(I,D) + kron(D,I);                % Kronecker sum since 2D
  Au = II + k*epsu*L; Av = II + k*epsv*L;          

% Time-stepping:
  for n = 1:1e6
    uold = u;
    u = Au*u + k*(-u.*v.^2 + F*(1-u));                  
    v = Av*v + k*(uold.*v.^2 - (c+F)*v);
    if mod(n,nplot)==0
      levs = 0:.1:1;
      uu = reshape(u,J,J);
      contourf(xx,yy,uu,levs)
      axis square, title(['t = ' num2str(k*n)], 'fontsize',12)
      axis on, drawnow
    end
  end