% solve the pde u_y +f(u,t,x) u_x = g(u,t,x) by method of characteristics

function pde_by_moc

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 
% set the x variable
% spatial grid
dx = 0.01;
x0 = [0:dx:1];
N = length(x0);
%initial data
u =x0.*(1-x0);
x = x0;
t_end = 3;
tstep = 0.25;
%specify the output points
tspan = [0:tstep:t_end];

s0=[x';u'];
 
[T,S] = ode23s(@deRHS,tspan, s0 );  

%plot the solution
% each row of S is a parametric representation of u, vs x
Tn = length(T);
N = length(x);
for j = 1:Tn
    plot(S(j,1:N),S(j,N+1:2*N),'linewidth',2)
    hold on
end



%the right hand side for ode simulation:
function s_prime=deRHS(t,s)

 N = length(s);
 x = s(1:N/2);
 u = s(N/2+1:N);
 
 % specify the MoC equations
 % u_t + Fx u_x = Fu
Fx = u;
Fu = zeros(N/2,1);
 
s_prime = [Fx; Fu];
