function traffic

% Define the regions over which the PDE is to be integrated
xmin = 0; xmax = 12; xint = .1;
tmin = 0; tmax = 12; tint = .1;
x = xmin: xint: xmax;
t = tmin: tint: tmax;

% Give the solver the various functions it needs in order to
% compute the PDE.
m = 0;
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);

u = sol(:,:,1);

surf(x,t,u)    
title('Unnecessary Shocks in Numerical Solution of PDE')
xlabel('Distance x (km)')
ylabel('Time t (min)')

r = 6;
c = 1;

% Plot selected density distributions at certain times
%for i=1:1:6,
%    subplot(6,1,i)
%    plot(x,u(i*10,:))
%    axis([0 12 0 10]);
%end

%for k=1:1:20
%    a(k) = max(integrate( fit(x', u(k,:)','cubicspline'),x,x(1)));
%end

%mina = min(a)
%maxa = max(a)

% Define the PDE in terms of the simple 1-D compressible gas model
function [c,f,s] = pdex1pde(x,t,u,DuDx)

MFS = .15;   % Define average velocity (mean free speed)

c = 1/(2*u);
f = u*MFS;
s = 0; 

% Create a type of density distribution for the t=0 initial
% conditions
function u0 = pdex1ic(x)
%u0 = 1/3 + 5/12*x;
u0 = 1 - 1/12*x;
%u0 = sin(pi*x/12);

% Give Neumann boundary conditions to the solver
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = ur;
qr = 0;