%Explicit numerical method for u_t + F(u) u_x = D u_xx

%%%%%%%%%% Numerical parameters
clear;  %  gets rid of old data
L = 150;  % domain size
eps = 0.2;   % "diffusion" constant - needed to smooth shocks
N = 1500; % number of spatial mesh points
dx = L/(N+1);
dt = 0.08*dx^2/eps^2;  %  timestep is limited by  1/2 (dx/eps)^2;
%setup grid
x = [dx: dx: L-dx];   % x is the column vector x_1, ..., x_N-1  % the ' denotes a matrix transpose
F = inline('.5 - 0.5*u','u','x');  %  define of F(u) in equation u_t + F(u) u_x = 0


%%%  initial data and dirichlet boundary conditions
%%%  Uncomment (remove the "%") to select;  be sure to comment the alternatives 

%%%  Step up initial data
u = 1*(x>20);  left = 0;  right = 1;

%%%  Step down initial data
%u = 1*(x<20);  left = 1;  right = 0;

%%%  "patch" initial data
%u = 2*rand(1,N).*(x>20 & x<30);  left = 0;  right = 0;

%%% no source term
source = inline('0','x','t');

%%% periodically varying source term at x=20
%source = inline('(.03+.03*sin(t/10))*exp(-(x-20).^2/5) ','x','t');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% step forward in time
time=0;  step = 0;
while (time < 2000),
   % these commands will plot the solution and pause
   % briefly before proceeding to the next time-step.
   % to make things go slower, increase the pause time
   if mod(step,100)==0, 
     plot(x,u);  axis([0 L 0 2]);   title(['time=' num2str(time)]); xlabel('x');  ylabel('u'); pause(.05);
   end;
     
   % do points next to boundary separately
   unext(1) =  u(1) + dt*( eps^2/dx^2*(-2*u(1) + u(2) + left )  - F(u(1),x(1))*(u(2)-left)/(2*dx)  );
   unext(N) = u(N) + dt*( eps^2/dx^2*(-2*u(N) + u(N-1) + right)   - F(u(N),x(N))*(right - u(N-1))/(2*dx) );
   %then do rest (note the "vectorization")
   unext(2:N-1) = u(2:N-1) + dt*(eps^2*(u(3:N)+u(1:N-2)-2*u(2:N-1))/dx^2 +... 
                    - F(u(2:N-1),x(2:N-1)).*(u(3:N)-u(1:N-2))/(2*dx) );
   
   u = unext + dt*source(x,time);          % now u is equal to the updated solution
   time = time + dt;  step = step+1;
end;


