% Example 13.7 Convection Equation
% Semi-discretization Method

clear all
close all
clc

% Physical parameters
c = 1;
L = 0.7;

% Mesh
dx = 0.01;
N = L/dx + 1;
x = 0:dx:L;

% time advancement
dt = 0.01;
tmax = 0.6;
t = 0:dt:tmax;
nmax = length(t);

u = zeros(nmax,N); % Allocate solution matrix

% initial condition
u(1,:) = exp(-200.*(x-0.25).^2);

% Define matrix
A = -1.*(c/2/dx).*(diag(zeros(1,N-1))+diag(-1.*ones(1,N-2),-1)+diag(ones(1,N-2),1));
A(N-1,N-2)=c/dx;
A(N-1,N-1)=-c/dx;

%Explicit Euler
for n = 1:nmax
    dudt = A*u(n,2:N)';
    u(n+1,2:N)=u(n,2:N)' + dt.*dudt;    
end

figure
%hold on
for n = 1:nmax
    plot(x,u(n,:),'-k')
    ylim([-4,4])
    xlabel('x')
    ylabel('u')
    title(['Explicit Euler solution of Convection Equation with dt = ',num2str(dt)])
    %saveas(gcf,[num2str(n),'.png'])
    pause(3*dt)
end

%tplots = [1, 401, 801, 1201, 1601];
tplots = [1, 13, 26, 39];
figure
hold on
for n = 1:length(tplots)
    plot(x,u(tplots(n),:),'LineWidth',2)     
end
ylim([-4,4])
xlabel('x')
ylabel('u')
title(['Explicit Euler solution of Convection Equation with dt = ',num2str(dt)])
legend(['t = ',num2str(t(1))],['t = ',num2str(t(13))],['t = ',num2str(t(26))],['t = ',num2str(t(39))])






