function characteristics( )
%CHARACTERISTICS
% Solves the equation
%   u_t + au_x = 0
% using the method of characteristics.

	% Set solver parameters
	nchar = 100;	% Number of characteristics
	T = 1;			% Maximum time
	nt = 100;		% Number of timesteps
	dom = [-3, 3];	% Initial spatial domain
	
	
	% Set the velocity function
	function a = velocity(u, x, t)
		a = 1;
% 		a = x.*(0<=x & x<1) + (1<=x);
% 		a = (x<0) + (1-x).*(0<=x & x<1);
% 		a = u;
	end


	% Set the initial data
	function u0 = initial(x)
		u0 = sin(pi*x);
% 		u0 = -atan(x);
	end


	% Set initial positions of characteristics, evaluate initial data and
	% set time parameters.
	xi = zeros(nt+1, nchar);
	xi(1,:) = linspace(dom(1), dom(2), nchar);
	u = initial(xi(1,:));
	dt = T / nt;
	t = linspace(0, T, nt+1);
	
	
	% Run simulation
	figure();
	for n = 2:nt+1
		% Compute positions of characteristics at next timestep
		xi(n,:) = xi(n-1,:) + dt*velocity(u, xi(n-1,:), t(n));
		
		% Plot solution
		subplot(1, 2, 1);
		plot(xi(1:n,:), t(1:n), 'b');
		ylim([0, T]);
		title('Characteristics');
		subplot(1, 2, 2);
		plot(xi(n,:), u, '-o');
		title(sprintf('u(t=%.2f)', t(n)));
		drawnow;
		
		% Halt execution after the first timestep
		if n == 2
			disp('Please press any key...');
			waitforbuttonpress();
			xlimits = xlim();
		end
	end
end