How to plot bifurcation diagram of the phased `1D chaotic map?

I plot bifurcation diagram of the map bellow but we get wrong bifurcation diagram not as authers plotted as picture bellow, please check and fix my code to get as the right bifurcation diagram , why i getted wrong bifurcation diagram?

Right bifurcation diagram:

The map is :

Di+1 =⎨16*G*Di*(0.5 − Di) 0 ≤ D < 0.25

16*G*(0.5 − Di)*(0.5 − G*(0.5 − Di)) 0.25 ≤ D < 0.5

16*G*(Di − 0.5)*(0.5 − G*(Di − 0.5)) 0.5 ≤ D < 0.75

6*G*(1 − Di)*(0.5 − G*(1 − Di))0.75 ≤ D ≤ 1

⎬

my code:

clear;

scale = 10000; % determines the level of rounding
maxpoints = 300; % determines maximum values to plot
N = 4000; % number of "r" values to simulate
a = 0; % starting value of "r"
b = 2; % final value of "r"... anything higher diverges.
rs = linspace(a,b,N); % vector of "r" values
M = 700; % number of iterations of  equation

% Loop through the "r" values
for j = 1:length(rs)
    
    G=rs(j); % get current "r"
    D=zeros(M,1); % allocate memory
    D(1) = 0.98; % initial condition (can be anything from 0 to 1)
    
    for i = 2:M, % iterate
        
    if(D(i-1) >=0 && D(i-1) <0.25)
        D(i)=16*G*D(i-1)*(0.5 - D(i-1));
    elseif (D(i-1) >=0.25 && D(i-1) <0.5)
        D(i)= 16*G*(0.5 - D(i-1))*(0.5 - G*(0.5- D(i-1)));
    elseif (D(i-1) >=0.5 && D(i-1) <0.75)
       D(i)= 16*G*(D(i-1) - 0.5)*(0.5 - G*(D(i-1) -0.5));
    elseif (D(i-1) >=0.75 && D(i-1) <=1)
        D(i)=16*G*(1 - D(i-1))*(0.5 - G*(1 - D(i-1)));
    else
        D(i)=0;
    end 
    end   

    % only save those unique, semi-stable values
    out{j} = unique(round(scale*D(end-maxpoints:end)));
end

% Rearrange cell array into a large n-by-2 vector for plotting
data = [];
for k = 1:length(rs)
    n = length(out{k});
    data = [data;  rs(k)*ones(n,1),out{k}];
end

% Plot the data
figure(2);clf
h=plot(data(:,1),data(:,2)/scale,'b.');
set(h,'markersize',1)
%axis tight
%set(gca,'units','normalized','position',[0 0 1 1])
%set(gcf,'color','white')
%axis off

my bifurcation diagram: