% code to solve Mackey-Glass delay differential equation
% using method of lines
function delay_de
global N beta dx

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 
formatSpecF = '%6.3f\n';
d = 5;  % time delay
betalist = [1/15;1/25];% decay time constant

for j = 1:length(betalist)
 beta = betalist(j);
N = 61;
dx = d/(N-1);
 

x = zeros(N,1); %initial data

tstep = 0.1;
t_end = 100;
tspan = [0:tstep:t_end];
 
s = zeros(N,1);
    
[T,S] = ode23s(@deRHS,tspan,s);
figure(j)
plot(T,S(:,1),'linewidth',2)
xlabel('time','fontsize',20)
ylabel('u(t)','fontsize',20)
 title(strcat('\beta = ',sprintf(formatSpecF,beta)),'fontsize',18)

end

function s_prime=deRHS(t,s)
global N beta dx
 
Fw = (s(1:N-1)-s(2:N))/dx;
Fs =  -beta*s(1) + 1/(1+s(N)^7);

s_prime = [Fs; Fw];

 




