Bifurcation code for delay Lorenz system

dx_1(t)/dt=&28*(x_2(t-\tau)-x_1(t)), dx_2(t)/dt=&10*x_1(t)-x_2(t)-x_1(t)x_3(t), dx_3(t)/dt=&-8/3*x_3(t-\tau)+x_1(t)x_2(t),

Here I want to draw a bifurcation diagram for x_1(t) vs. \tau. \tau is taken as bifurcation parameter. My code is given below. Please check and correct it. I am not getting the bifurcation diagram.

MATLAB CODE:

function rossler;

% 3-variable Rossler model - chaos % Didier Gonze % 8/7/2008

clc;

%%%% Number of variable and initial conditions:

% nbvar=3; % xini=ones(1,nbvar)/10; % % %%%% Time parameters: % trans=100; tend=500; tstep=0.01; %

history = [0.1; 0.1; 0.1]; %%%% Range (for bifurcation diagram as a function of b):

b=0.2; % (default value for chaos) Z=b; bmin=0.2; bmax=0.3; bint=0.05;

brange=[bmin bint bmax];

%%%% Task:

%integration(xini,trans,tend,tstep,b); bifurcation(brange,history,trans,tend,tstep);

%==================================================================== % Integration %====================================================================

% function output=integration(x0,trans,tend,tstep,b); % % [t,x] = run(x0,trans,tend,tstep,b); % % % set(figure(1),'Position', [400 400 500 300]); % clf; % % plot(t,x(:,1:3)); % xlabel('Time','fontsize',18); % ylabel('Variables','fontsize',18); % xlim([0 tend]); % legend('X','Y','Z'); % % % set(figure(2),'Position', [400 400 500 300]); % clf; % % plot3(x(:,1),x(:,2),x(:,3)); % xlabel('X','fontsize',18); % ylabel('Y','fontsize',18); % zlabel('Z','fontsize',18); % box on; %

%==================================================================== % Bifurcation %====================================================================

function output=bifurcation(range,history,trans,tend,tstep)

D=[]; % data (bifurcation diagram)

for Z=range(1):range(2):range(3)

    fprintf('Z=%g...\n',Z);

    sol = run(Z,history,trans,tend,tstep);

    for i=2:length(sol.y(:,2))-1
        if((sol.y(i,2)>sol.y(i-1,2))&&(sol.y(i,2)>sol.y(i+1,2)))
            D=[D; Z sol.y(i,2)];
        end
    end

end

figure(3) plot(D(:,1),D(:,2),'ro','MarkerEdgeColor','b','MarkerFaceColor','b','MarkerSize',1.5)

xlabel('b','fontsize',18); ylabel('max(X)','fontsize',18);

% =================================================================== % Run % ===================================================================

function sol=run(Z,history,trans,tend,tstep)

ttrans = [0:tstep:trans]; tspan = [0:tstep:tend];

option = []; %option = odeset('RelTol', 1e-5); %option=odeset('OutputS',[1:3],'OutputF','odeplot');

if trans > 0 sol = dde23(@dxdt,Z,history,ttrans); x0=sol.y(end,:); end

sol = dde23(@dxdt,Z,history,tspan);

% =================================================================== % dxdt % ===================================================================

function y = dxdt(t,x,Z)

%%% parameters y=zeros(3,1); xlag=Z(:,1); %delay

%%% equations

y = [ 10*(xlag(2)-x(1)); % dx/dt x(1)*(28-x(3))-x(2); % dy/dt x(1)*x(2)-8/3*xlag(3); % dz/dt ] ;