% math4090 2010 exercise problem
% discrete traffic dynamics - delayed differential equation
%
clear all
%
v=200e3/3600;
rho_max=40e-3;
L=1/rho_max;
rho_o=rho_max*exp(-1);
d=1/rho_o-L;
N=5;
T=50;
NT=501;
dt=T/(NT-1);
t=[0:dt:T];
tau=1.1;
Ntau=int16(tau/dt);
z=zeros(NT,N);
t0=1;
k=0.2;
z(:,1)=(v*exp(1)*k*t0)*((t+t0).*exp(-t/t0)-t0);
%
for j=2:N
    Nj=1+(j-1)*Ntau;
    for i=Nj:NT-1
        rhs=v*log(1+rho_max*exp(-1)*(z(i-Ntau,j-1)-z(i-Ntau,j)));
        z(i+1,j)=z(i,j)+dt*rhs;
    end
end
%
y=zeros(NT,N);
for j=2:N
    y(:,j)=-(j-1)*(d+L);
end
x=y+z;
%
figure(1);
plot(t,x(:,:));
xlabel('t');
ylabel('z_i(t)');