% Example 13.2 Convection-Diffusion equation
% Using FTCS Method

clear all
close all
clc

% Physical parameters
L = 1.0;
alpha = 0.010;
u = 0.10;
tmax = 50;
bc1 = 0;
bc2 = 100;

% Mesh
dx = 0.1;
dt = 0.5;
c = u*dt/dx;
d = alpha*dt/dx^2;

imax = L/dx+1;
nmax = tmax/dt+1;
f = zeros(nmax,imax); % allocate solution matrix
x = [0:dx:L];
t = [0:dt:tmax];

% initial condition
f(1,:) = 100.*x./L;

% boundary conditions
f(:,1) = bc1;
f(:,imax) = bc2;

% FTCS Method
for n = 1:nmax-1
    for i = 2:imax-1
        f(n+1,i) = f(n,i) - (c/2)*(f(n,i+1)-f(n,i-1)) + d*(f(n,i+1)-2*f(n,i) + f(n,i-1));
    end
end

figure
plot(x,f(1,:),'-','LineWidth',2)
xlabel("x (cm)")
ylabel("T ^oC")
title("Initial temperature profile")

% Plot several lines of the solution
tsol = [0.5, 2.5, 5.0, 10, 50];
nsol = tsol./dt+1;
marks = ['s','^','o','>','*'];

% Analytical Solution
xexact = [0:0.00001:L];
for j = 1:length(nsol)
    for i = 1:length(xexact)
        fexact(j,i) = convdiff(u,alpha,L,xexact(i),tsol(j));
    end
end

figure
hold on
plot(x,f(1,:),'-k','LineWidth',2)
for j = 1:length(nsol)
    plot(x,f(nsol(j),:),marks(j))
end
for j = 1:length(nsol)
    plot(xexact,fexact(j,:),'-')
end
ylim([0,100])
xlabel('x')
ylabel('T ^oC')
title("Convection-Diffusion Equation - FTCS Method")
legend('t=0','t=0.5','t=2.5','t=5','t=10','t=50','t=0','t=0.5','t=5','t=10','t=50','Location','northwest')


figure
ylim([0,100])
xlabel('x')
ylabel('T ^oC')
title("Convection-Diffusion Equation - FTCS Method")
hold on
for j = 1:nmax/3
    plot(x,f(j,:),'-','MarkerEdgeColor',[0.8500, 0.3250, 0.0980])
    pause(dt);
end

figure
surf(x,t,f)
xlabel('x')
ylabel('t')
zlabel('T ^oC')
