function [x,y,u] = RDzeroflux_spots(N,a,b,d,gamma)
% defining the mesh
N = 30;
h = 1/N; % step size in x and y
x = h*(0:N); % x coordinates of grid
y = h*(0:N); % y coordinates of grid
[xx,yy] = meshgrid(x,y); % 2D x and y mesh coordinates
dt = .01*h^2; % time step - usually small
% parameter values
a = .05;
%b = 1.4;
b = 1;
%d = 14.85;
d = 20;
%gamma = 400;
gamma = 600;
% Initial data at t=0:
us = (a+b); % u steady state
vs = b/us^2; % v steady state
u = (a+b)*ones(size(xx)); % u steady state
v = (b./u.^2); % v steady state
u = u + .1*randn(size(xx)); % add small perturbations about the steady state
v = v + .1*randn(size(xx)); % add small perturbations about the steady state
% Time-stepping:
t = 0;
tmax = 10;
nsteps = round(tmax/dt); % number of time steps
wb = waitbar(0,'RD Spots');
nsteps = min(50000,nsteps);
for n = 1:nsteps % main time-stepping loop
    waitbar(n/nsteps);
t = t+dt;
uE = u(:,[2:N+1 N]);
uW = u(:,[2 1:N]);
uN = u([2 1:N],:);
uS = u([2:N+1 N],:);
vE = v(:,[2:N+1 N]);
vW = v(:,[2 1:N]);
vN = v([2 1:N],:);
vS = v([2:N+1 N],:);
% finite difference formula
u2v = u.^2.*v;
u = u + gamma*dt*(a - u + u2v) + dt*(uE+uW+uN+uS-4*u)/h^2;
v = v + gamma*dt*(b - u2v) + d*dt*(vE+vW+vN+vS-4*v)/h^2;
end
close(wb);
% plot solution:
colormap('gray')
surf(x,y,u)
title(['u at t = ' num2str(t)],'fontsize',16)
zlim([0 4])
shg

