function [X,v] = quasilin(c,f,u0,xval,tval)
% [X,v] = quasilin(c,f,u0,xval,tval)
% solve quasilinear 1st order PDE 
%   u_t + c(x,t,u)*u_x = f(x,t,u),    u(0)=u0(x)
% using method of characteristics (solving ODEs with ode45)
% 
% inputs:
%   c:    coefficient fct c, define by c = @(x,t,u) ...
%   f:    right hand side fct f, define by f = @(x,t,u) ...
%   u0:   initial condition fct u0, define by u0 = @(x) ...
%   xval: values for x0, define by xval = xmin:step:xmax
%   tval: values for t,  define by tval = tmin:step:tmax
% outputs: 
%   X: x-values of characteristics
%   v: u-values of solution
% X,v are arrays where row index corresponds to x0-values xval
%                   column index corresponds to t-values tval
% at time tval(j) the x-values are X(:,j), the u-values are v(:,j)
% to plot the solution u vs. x at time tval(3) use  plot(X(:,3),v(:,3))
%
% Example: Burgers equation u_t + u*u_x = 0, u(x,0)=atan(x)
%   c = @(x,t,u) u;       % c(x,t,u)=u
%   f = @(x,t,u) 0;       % f(x,t,u)=0
%   u0 = @(x) atan(x);    % initial condition u0(x)=atan(x)
%  [X,v] = quasilin(c,f,u0,-5:.1:5,0:.1:3);

% This uses colorcurves.m 
% Download colorcurves.m from the course web page!

Nx = length(xval);
Nt = length(tval);
X = zeros(Nx,Nt);
v = zeros(Nx,Nt);
%  need to solve ODE system  X'=c(X,t,v), v'=F(X,t,v)
%  let w = (X,v), i.e., w(1)=X, w(2)=v
%  then we can write ODE system as w'=G(t,w)
G = @(t,w) [ c(w(1),t,w(2)) ; f(w(1),t,w(2)) ]; % define fct G
for i=1:Nx
  x0 = xval(i); 
  [ts,ws] = ode45(G,tval,[x0;u0(x0)]); % solve system of ODEs w'=G(t,w)
  X(i,:) = ws(:,1);   % 1st column of ws are values of X, store as row
  v(i,:) = ws(:,2);   % 2nd column of ws are values of v, store as row
end
% each row of array X contains solution for one value of x0
% each row of array v contains solution for one value of x0

% plot characteristics in (x,t) plane
figure(1)
plot(X',tval,'k')     % plot rows of X vs tval
xlabel('x'); ylabel('t')
title('characteristics')

% plot solutions u(x,t_j) vs x, for all times
figure(2)
colorcurves(X,tval,v)  % plot columns of v vs. columns of X
xlabel('x'); ylabel('t'); zlabel('u')
title('Rotate this graph with the mouse')




