% to solve u_t + J_x= 0, with MoL and upwinding (written in conservation
% form, J = v*u

%this is for an exercise with burgers equation u_t + (u^2)_x=0
% J=u^2
function pde_MoL_upwind
global dx N x

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 

% spatial grid
dx = 0.01;
x = [-5:dx:5];
N = length(x);
%initial data
u =1./(1+x.^2);


% try method of lines (MoL) with upwinding
t_end = 2;  % 
tstep = 0.2;
%specify the output points
tspan = [0:tstep:t_end];

u0 = u'; % must be a column vector
%solve the  MOL differential equations
[T,S] = ode23s(@deRHS,tspan, u0 );  

 
figure(1)
 plot(x,S(1,:),'--','linewidth',2)
  hold on
for j = 2:length(T)
    plot(x,S(j,:),'linewidth',2)
  hold on
end
xlabel('x','fontsize',20)
ylabel('u(x,t)','fontsize',20)
set(gca,'Xtick',-5:2.5:5)
hold off

% the solution by method of characteristics:
x0 = [-5:dx:5];
u0 =1./(1+x0.^2);

figure(2)
    plot(x0,u0,'--','linewidth',2)
    hold on
for j = 2:length(T)
    t = T(j);
    x = 2*t.*u0+x0;
    figure(2)
    plot(x,u0,'linewidth',2)
    hold on
end
xlabel('x','fontsize',20)
ylabel('u(x,t)','fontsize',20)
set(gca,'Xtick',-5:2.5:5)
axis([-5 5 0 1.2])
hold off

% plot the cusp
[x,t]=meshgrid(0:.05:5,0:0.05:5);
R= 8*t.*x.^3 - 4*x.^4 - 108*t.^2 + 72*t.*x - 8*x.^2 - 4;
figure(3)
contour(x,t,R,[0 0],'linewidth',2)
xlabel('x','fontsize',20)
ylabel('t','fontsize',20)
text(1,3,'I','fontsize',20)
text(4,3,'II','fontsize',20)


%the right hand side for ode simulation:
function s_prime=deRHS(t,u)
global dx N x

% by definition, J = v*u, where v = v(x,t,u) in general
% first evaluate v_{j-1/2}
 
xjmh = [x-dx/2,x(N)+dx/2]'; % (x_{j-1/2})
ujmh = ([u;0]+[0;u])/2;  %u_{j-1/2} The zero fill is a ghost point outside the grid

%use xjmh and ujmh to evaluate v_{j-1/2} = v(xjmh,t,ujmh) 
vjmh =   ujmh;  %this is for this exercise, J = u^2

 %this is upwinding:
Jmh = vjmh.*((vjmh>0).*[0;u]+(vjmh<0).*[u;0]);

Fu = (Jmh(1:end-1)-Jmh(2:end))/dx;  %finite difference in x
 
 
s_prime =  Fu;