function [x,y,u] = RDzeroflux_spots2(N,a,b,d,gamma,nei,d_max)
% defining the mesh
N = 30;
[num_pts,num_pts] = size(nei);
num_side = sqrt(num_pts);
nnei = nei<d_max&nei>0;
h = 1/(num_side-1); % step size in x and y
x = h*(1:num_side); % x coordinates of grid
y = h*(1:num_side); % 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(5000,nsteps);
for n = 1:nsteps % main time-stepping loop
%    [n nsteps]
    waitbar(n/nsteps);
    t = t+dt;
    udiff = zeros(1,num_pts);
    vdiff = zeros(1,num_pts);
    for p = 1:num_pts
        row = rem(p-1,num_side) + 1;
        col = floor((p-1)/num_side) + 1;
        udiff(p) = sum(nnei(p,:).*u(:)') - sum(nnei(p,:).*u(row,col));
        vdiff(p) = sum(nnei(p,:).*v(:)') - sum(nnei(p,:).*v(row,col));
    end
    udiff = reshape(udiff,num_side,num_side);
    vdiff = reshape(vdiff,num_side,num_side);
    % finite difference formula
    u2v = u.^2.*v;
    u = u + gamma*dt*(a - u + u2v) + dt*(udiff)/h^2;
    v = v + gamma*dt*(b - u2v) + d*dt*(vdiff)/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
