% Example 13.7 Convection Equation
% Semi-discretization Method

clear all
close all
clc

% Physical parameters
c = 1;
L = 1;

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

% time advancement
dt = 0.01;
tmax = 1.1;
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;

%RK4
for n = 1:nmax
    k1 = A*u(n,2:N)';
    k2 = A*u(n,2:N)' + dt/2.*A*k1;
    k3 = A*u(n,2:N)' + dt/2.*A*k2;
    k4 = A*u(n,2:N)' + dt.*A*k3;
    
    u(n+1,2:N)=u(n,2:N)' + dt.*(k1+2.*k2+2.*k3+k4)./6;   
end


figure
%hold on
for n = 1:nmax
    plot(x,u(n,:),'-k')
    ylim([-0.1,1.5])
    yticks([0,0.25,0.5,0.75,1,1.25,1.5])
    xticks([0,0.25,0.5,0.75,1])
    xlabel('x')
    ylabel('u')
    title(['RK4 solution of Convection Equation with dt = ',num2str(dt)])
    %saveas(gcf,[num2str(n),'.png'])
    pause(dt)
end

tplots = [1, 26, 51, 76];
figure
hold on
for n = 1:length(tplots)
    plot(x,u(tplots(n),:),'LineWidth',2)     
end
ylim([-0.1,1.5])
yticks([0,0.25,0.5,0.75,1,1.25,1.5])
xticks([0,0.25,0.5,0.75,1])
xlabel('x')
ylabel('u')
title(['RK4 solution of Convection Equation with dt = ',num2str(dt)])
legend(['t = ',num2str(t(1))],['t = ',num2str(t(26))],['t = ',num2str(t(51))],['t = ',num2str(t(76))])






