function pdex6
  
% In the form expected by PDEPE, the PDE is
%
%    [1] .*  D_ [p] =  D_ [   0   ] +  [ -v*Dp/Dx ]
%            Dt        Dx
%    ---        ---        -------------      ----------------
%     c          p         f(x,t,p,Dp/Dx)      s(x,t,p,Dp/Dx)
%
% Here v = dq/dp = p*du/dp + u in the traffic flow model, with p
% representing the density, q representing the flux, and u the velocity
% which is either constant or linear as determined below in pdex6pde.
%
% MATLAB requires a finite domain, thus 0<=x<=b, and 0<=t<=T.
%
% The initial condition is given by 
%
%   p(x,0) = p0(x)
%
% where f(x) is defined as an internal function with examples below.
%
% The left boundary condition is p(0,t) = 0;
%
%      [p]    +     [0] .* [       Dp/Dx          ] = [0]
%
%      ---          ---    ------------------------   ---
%    p(0,t,p)      q(0,t)        f(0,t,p,Dp/Dx)        0
%
%   The right bc is p(L,t) = 0:
%
%      [p]    +     [0] .* [       Dp/Dx          ] = [0]
%
%      ---          ---    ------------------------   ---
%    p(L,t,p)      q(L,t)        f(L,t,p,Dp/Dx)        0
%
  
% Problem parameters, shared with nested functions.
  b=10;
  T=10;
  umax=.5;
  pmax=10;
  
  m = 0;
  n = 101; % increase this number for better numerical resolution
  x = linspace(0,b,n);
  t = linspace(0,T,21);
  p = pdepe(m,@pdex6pde,@pdex6ic,@pdex6bc,x,t);
  figure;
  surf(x,t,p)
  title(['Numerical solution computed with ', num2str(n), ' mesh points.'])
  xlabel('Distance x')
  ylabel('Time t')
  zlabel('Density \rho')
  rotate3d on
  
  figure
  plot(x,p(1,:),'-o',x,p(end,:),'--*')
  title(['Initial and final density profiles'])
  xlabel('Distance x')
  ylabel('Density \rho')
  legend('Initial','Final')
  
% -----------------------------------------------------------------------
% Nested functions -- parameters are provided by the outer function.
%
  
function [c,f,s] = pdex6pde(x,t,p,DpDx)
% Evaluate the differential equations components.
  c = 1;                 
  f = 0*DpDx;            
  u = umax; dudp = 0; % constant velocity-density relationship
  %u = umax*(1-p/pmax); dudp = -umax/pmax; % linear velocity-density
  s = -(p*dudp + u)*DpDx;       
  end
% -----------------------------------------------------------------------

function p0 = pdex6ic(x)
  p0 = exp(-(x-b/4)^2);
  %p0 = (x<4).*x;
  %p0 = x>1 & x<4; or try %p0 = exp(-((x-b/4)/1.5)^100);
  end
% -----------------------------------------------------------------------

function [pl,ql,pr,qr] = pdex6bc(xl,pl,xr,pr,t)
% Evaluate the boundary conditions components.
  pl = pl;                 
  ql = 0;                 
  pr = pr;                 
  qr = 0;                 
  end 
% -----------------------------------------------------------------------

end  % pdex6