function [x,y,u] = RDzeroflux_spots_err(N,a,b,d,gamma,s_scale)
% 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

% Displace points
for xi = 1:N+1
    for yi = 1:N+1
        xv = xx(xi,yi);
        yv = yy(xi,yi);
        x_del = 0.5 - rand;
        y_del = 0.5 - rand;
        xx(xi,yi) = xv + s_scale*x_del;
        yy(xi,yi) = yv + s_scale*y_del;
    end
end

% Distances between neighbors
dE = zeros(N+1,N+1);
for r = 1:N+1
    for c = 1:N
        dE(r,c) = norm([xx(r,c),yy(r,c)]-[xx(r,c+1),yy(r,c+1)]);
    end
    dE(r,N+1) = norm([xx(r,N+1)-1,yy(r,N+1)]-[xx(r,1),yy(r,1)])+h;
end
dW = dE(:,[N+1 1:N]);
dS = zeros(N+1,N+1);
for c = 1:N+1
    for r = 1:N
        dS(r,c) = norm([xx(r,c),yy(r,c)]-[xx(r+1,c),yy(r+1,c)]);
    end
    dS(N+1,c) = norm([xx(1,c),yy(1,c)]-[xx(N+1,c),yy(N+1,c)-1])+h;
end
dN = dS([N+1 1:N],:);

if s_scale==0
    dE = ones(N+1,N+1)*h^2;
    dW = dE;
    dS = dE;
    dN = dE;
end

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');
num2do = min(nsteps,50000);
for n = 1:num2do % main time-stepping loop
    waitbar(n/num2do);
    t = t+dt;
    uE = u(:,[2:N+1 N]);
%    uW = u(:,[2 1:N]);
    uW = u(:,[N+1 1:N]);
%    uN = u([2 1:N],:);
    uN = u([N+1 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;
    Lap_u1 = uE + uW + uN + uS - 4*u;
    Lap_u = ((uW - u)./dW - (u - uE)./dE) + ((uN - u)./dN - (u - uS)./dS);
    Lap_v1 = vE + vW + vN + vS - 4*v;
    Lap_v = ((vW - v)./dW - (v - vE)./dE) + ((vN - v)./dN - (v - vS)./dS);
%    u = u + gamma*dt*(a - u + u2v) + dt*(Lap_u)/h^2;
%    v = v + gamma*dt*(b - u2v) + d*dt*(Lap_v)/h^2;
    u = u + gamma*dt*(a - u + u2v) + dt*(Lap_u);
    v = v + gamma*dt*(b - u2v) + d*dt*(Lap_v);
    if rem(n,1000)==0
        tch = 0;
    end
end
close(wb);
% plot solution:
colormap('gray')
surf(x,y,u)
title(['u at t = ' num2str(t)],'fontsize',16)
zlim([0 4])
shg
