AE5139HW5

% AE5139 - Intro to CFD
% Homework 5
% Due: 19 September 2019
% Christopher Bates
% Pastebin Link:

close all
clear all
clc
trep = [0.05,0.1,0.2,0.5,5];                 % reported time steps

t0               = 0;
tf               = 5;
numtimesteps     = 20;                       % Short program, so 100 shouldnt take long
numintervals     = 20;                       % number of x intervals
L                = 1;                        % Given length
dx               = L/numintervals;
xloc             = 0:dx:L;                   % Generating mesh
dy               = 0;                        % infinite slab, 1-D HT

dt       = 0.001;

% Ensuring dt is small enough for stability
if dt > 0.5*(dx^2 + dy^2)
  dt = 0.5*(dx^2 + dy^2);
end

lambda  = dt/dx^2;
t       = dt:dt:tf;

% Initialize mesh
u = zeros(length(t),length(xloc));

% Set starting conditions on mesh
u(1,:) = 250*sin(pi.*xloc);                     % K, Given starting temp function
u(1,1)              = 350;                      % K, Boundary condition at x = 0
u(1,numintervals+1) = 420;                      % K, Boundary condition at x = L0 = xmax

for i = 2:length(t)
    u(i,1)              = 350;
    u(i,numintervals+1) = 420;
    for j = 2:numintervals
        u(i,j) = lambda*u(i-1,j-1) + (1-2*lambda)*(u(i-1,j))...
            + lambda*u(i-1,j+1);
    end
end

figure
hold on
for k = trep
    txt = ['t = ',num2str(k)];
    plot(xloc,u(k/dt,:),'DisplayName',txt)
end
grid on
xlabel('Nondimensional Length','Interpreter','latex')
ylabel('Temperature (K)','Interpreter','latex')
title('Temperature Vs. Non-dimensional Length (Explicit Method)','Interpreter','latex')
legend show
legend('Location','best')