clear;
% timestep
DeltaT=0.1;
% diffusion coefficients
DX=0.2;
DY=0.02;
% reaction constants
k1=1;
k2=1;
k3=1;
k4=1;
% initial concentrations
X0=1;
Y0=1;
% constant concentrations
A=1;
B=3;
% sigma of the random initial distortion
initRandomSigma=2;
% width/height of the simulation grid
width=100;
% number of iterations before display
iterationsPre=3000;
% number of visible iterations after pre
iterations=300;
iterationsDT=6;
% number of visible iterations before pre
iterationsStart=300;
% number of visible iterations before pre
iterationsInit=300;
% number of invisible iterations before iterationsStart
iterationsStartPre=600;
% number of visible iterations after pre
iterationsPre1=6000;
% number of visible iterations after pre
iterations1=300;
 
 
% laplace operator
laplace=0.25*[0 1 0; 1 -4 1; 0 1 0];
 
% initialize reaction variables
X(1:width,1:width)=X0+(rand(width,width)-0.5)*2*initRandomSigma;
Y(1:width,1:width)=Y0+(rand(width,width)-0.5)*2*initRandomSigma;
X(X<0)=0;
Y(Y<0)=0;
X(X>4.5)=4.5;
Y(Y>4.5)=4.5;
 
Xi=X;
Yi=Y;
 
% colormap
rbmap(1:100,1:3)=0;
rbmap(:,1)=((length(rbmap)-1)-(0:(length(rbmap)-1)))/(length(rbmap)-1);
rbmap(:,3)=(0:(length(rbmap)-1))/(length(rbmap)-1);
 
 
figure(1)
s=get(gcf, 'Position');
s(3)=250;
s(4)=200;
set(gcf, 'Position', s);
mf=1;
for it=1:(iterationsInit+iterationsStartPre+iterationsStart+iterationsPre+iterations+iterationsPre1+iterations1)
    lapX=conv2(X, laplace, 'same');
    lapY=conv2(Y, laplace, 'same');
    vX=k1*A+k2*X.^2.*Y-k3*B*X-k4*X+DX.*lapX;
    vY=k3*B*X-k2*X.^2.*Y+DY.*lapY;
 
    % leave border as distrotion centers
    vX(1:width,[1 width])=0;
    vX([1 width],1:width)=0;
    vY(1:width,[1 width])=0;
    vY([1 width],1:width)=0;
 
    if ((it<=iterationsInit) ...