I want to find basin of attraction for delayed logistic map like this:

x(n+1) = A*x*(1-x(n)-B*x(n-1)

I wrote bellow (MATLAB) code, but it just show three basin of attraction for given parameter value.

%**************************************************************************************

clc ; clear all, close all

a = 3.815; % parameter

b = 0.19; % parameter

n = 1500; % Number of iterations

% Initial conditions for the Attractor

X(1) = 0.037;

Y(1) = 0.037;

% Apply the iterations

for ii = 1:n

X(ii+1) =a*X(ii)*(1 -X(ii)) -b* Y(ii);

Y(ii+1)= X(ii);

end

% Initial conditions for the Attractor2

X1(1) = 0.8;

Y1(1) = 0.8;

% Apply the iterations

for ii = 1:n

X1(ii+1) =a*X1(ii)*(1 -X1(ii)) -b* Y1(ii);

Y1(ii+1)= X1(ii);

end

% Initial conditions for the basin of attraction

xx = linspace(0,1,100) ;

yy = linspace(0,1,100) ;

% Tolerance for considering convergence

tolerance = 1e-2;

% Initialize the required matrices

Attract1 = [] ; Divergence = [] ; Attract2 = [];c=[];

tic

for i = 1:length(xx)

for j = 1:length(yy)

X0 = [xx(i);yy(j)];

x = xx(i);

y = yy(j);

% Apply the iterations

for k = 1:n

x_next =a *x*(1 - x)-b*y;

y = x;

x = x_next;

end

% Classify the point based on its final location

if norm(min (abs([X' Y']-[x*ones(size(X))' y*ones(size(Y))']))) < tolerance

Attract1 = [X0 Attract1] ;

elseif norm(min (abs([X1' Y1']-[x*ones(size(X1))' y*ones(size(Y1))']))) < tolerance

Attract2 = [X0 Attract2] ;

else % if not close to the attractor

Divergence = [X0 Divergence] ;

end

end

end

toc

%Initialize figure

figure

set(gcf,'color','w')

hold on

plot(Attract1(1,:),Attract1(2,:),'.','color','r') ;

plot(Attract2(1,:),Attract2(2,:),'.','color','g') ;

plot(Divergence(1,:),Divergence(2,:),'.','color','b') ;

% plot(X(100:end),Y(100:end),'r.','MarkerSize',4)

title('Basin of attraction')