
%% Grid for periodic boundary conditions
Nx = 60;  hx = 1/Nx;  x1 = (hx:hx:1)';
Ny = 60;  hy = 1/Ny;  y1 = (hy:hy:1)';
[xx,yy] = meshgrid(x1,y1);
x = xx(:); y = yy(:);

e = ones(Nx,1);
Dx = spdiags([e  -2*e  e], [-1 0 1], Nx, Nx);
Dx(1,end) = 1;  Dx(end,1) = 1;  % periodic BCs
Dx = 1/hx^2 * Dx;
Ix = speye(size(Dx));

e = ones(Ny,1);
Dy = spdiags([e  -2*e  e], [-1 0 1], Ny, Ny);
Dy(1,end) = 1;  Dy(end,1) = 1;  % periodic BCs
Dy = 1/hy^2 * Dy;
Iy = speye(size(Dy));
Dxx = kron(Dx, Iy);
Dyy = kron(Ix, Dy);
L = Dxx + Dyy;
uinitf = @(x,y) min(1,10*sqrt((x-0.2).^2+(y-0.2).^2));
vinitf = @(x,y) max(0,1-10*sqrt((x-0.3).^2+2*(y-0.3).^2));
%u = uinitf(x,y);
%v = vinitf(x,y);
k = 2.5*hx;

Tf = 3000;
numsteps = ceil(Tf / k);
k = Tf / numsteps;
coords1 = find(x1==0.75);
coords2 = find(y1==0.75);
format long
%
%%
% Case 1, t=500 and t=3000
c = 0.065;
F = 0.06;
u = uinitf(x,y);
v = vinitf(x,y);
coeff1 = 6e-5;
coeff2 = 2e-5;
for n = 1:numsteps
  u = u + k*(coeff1*L*u-u.*v.^2+F*(1-u));
  v = v + k*(coeff2*L*v+u.*v.^2-(c+F)*v);
  if mod(n,50)==0
    uu = reshape(u,Nx,Ny);
    pcolor(xx,yy,uu);
    %waitforbuttonpress;
    drawnow;
  end
end
%}

%{
%%
% Case 2, t = 500 and t=3000
figure
c = 0.065;
F = 0.03;
u = uinitf(x,y);
v = vinitf(x,y);
for n = 1:numsteps
  u = u + k*(coeff1*L*u-u.*v.^2+F*(1-u));
  v = v + k*(coeff2*L*v+u.*v.^2-(c+F)*v);
  if mod(n,50)==0
    uu = reshape(u,Nx,Ny);
    pcolor(xx,yy,uu);
    drawnow;
  end
end
%}