sample

clc;
clear;

%% Plate Dimensions
Ly = 3.5; %Height (in)
Lx = 0.3125; %Width (in)

%% Spatial Discretization of Elements
dx = 0.0521; %(in)
dy = 0.140; %(in)
dt = 0.08; %(s)
n_i = 26; %Number of nodes along the i-axis
n_j = 7; %Number of nodes along the j-axis

%% Material Properties and Constants
K = 1.9907e-4; %Conductivity of Steel (Btu/hr-in-R)
rho = .2853; %Density (lbm/in^3)
c_p = 0.114; %Specific Heat (Btu/lbm-R)
h = 3.858e-5; %Heat Transfer Coefficient (Btu/hr-in^2-F)
a = 6.12e-3; %Thermal Diffusivity (in^2/s)

%% Inital Conditions
T_inf = 60; %Air temperature (F)
Ts = T_inf; %Surface temperature (F)
q_dot = 0.020; %Heat Flux (Btu/s-in^2)
T(1:n_i,1:n_j,1)= Ts; %All nodes along the right wall are at Ts
t = 0; %Initial time (s)
t_f = 8; %Final stopping time (s)

%% Temperature Profile Generation

k = 1; %Beginning iteration at the first time step

while t < t_f
    for i = 1:n_i
        for j = 1:n_j

% General Recursive Formula for Interior Nodes
    if (i>1 && i<n_i) && (j>1 && j<n_j)
          T(i,j,k+1)= T(i,j,k)*(1 - (a*dt)/(dx^2) -(a*dt)/(dx^2) - (a*dt)/(dy^2) - (a*dt)/(dy^2)) + T(i,j-1,k)*(a*dt)/(dx^2) + T(i,j+1,k)*(a*dt)/(dx^2) + T(i+1,j,k)*(a*dt)/(dy^2)+ T(i-1,j,k)*(a*dt)/(dy^2);

% Top nodes (52,78,104,130,156)
        elseif i == n_i && (j>1 && j<n_j)
                T(i,j,k+1)= T(i,j,k)*(1 - (2*h*dt)/(rho*c_p*dy) -(a*dt)/(dx^2) - (a*dt)/(dx^2) - (2*a*dt)/(dy^2)) + T_inf*(2*h*dt)/(rho*c_p*dy) + T(i,j-1,k)*(a*dt)/(dx^2) + T(i,j+1,k)*(a*dt)/(dx^2) + T(i-1,j,k)*(2*a*dt)/(dy^2);

% Bottom nodes (27,53,79,105,131)
        elseif i == 1 && (j>1 && j<n_j)
                T(i,j,k+1)= T(i,j,k)*(1 - (a*dt)/(dx^2) -(a*dt)/(dx^2) - (2*a*dt)/(dy^2)) + T(i,j-1,k)*(a*dt)/(dx^2) + T(i,j+1,k)*(a*dt)/(dx^2) + T(i+1,j,k)*(2*a*dt)/(dy^2);

% Right nodes (158,159,160, ... 181)
        elseif j == n_j && (i>1 && i<n_i)
                T(i,j,k+1)= T(i,j,k)*(1 - (2*a*dt)/(dx^2) -(2*h*dt)/(rho*c_p*dx) - (a*dt)/(dy^2) - (a*dt)/(dy^2)) + T_inf*(2*h*dt)/(rho*c_p*dx) + T(i,j-1,k)*(2*a*dt)/(dx^2) + T(i+1,j,k)*(a*dt)/(dy^2) + T(i-1,j,k)*(a*dt)/(dy^2) +((2*q_dot*dt)/(rho*c_p*dx));

% Left nodes (2,3,4, ... 25)
        elseif j == 1 && (i>1 && i<n_i)
                T(i,j,k+1)= T(i,j,k)*(1 - (2*h*dt)/(rho*c_p*dx) -(2*a*dt)/(dx^2) - (a*dt)/(dy^2) - (a*dt)/(dy^2)) + T_inf*(2*h*dt)/(rho*c_p*dx) + T(i,j+1,k)*(2*a*dt)/(dx^2) + T(i+1,j,k)*(a*dt)/(dy^2) + T(i-1,j,k)*(a*dt)/(dy^2);

% Upper Left Node (1)
        elseif i==1 && j==1
          T(i,j,k+1)= T(i,j,k)*(1 - (2*h*dt)/(rho*c_p*dx) - (2*a*dt)/(dx^2) - (2*a*dt)/(dy^2)) + T_inf*((2*h*dt)/(rho*c_p*dx)) + T(i,j+1,k)*(2*a*dt)/(dx^2) + T(i+1,j,k)*(2*a*dt)/(dy^2);

% Lower Left Node (26)
        elseif j==1 && i==n_i
          T(i,j,k+1)= T(i,j,k)*(1 - (2*h*dt)/(rho*c_p*dx) -(2*h*dt)/(rho*c_p*dy) - (2*a*dt)/(dx^2) - (2*a*dt)/(dy^2)) + T_inf*((2*h*dt)/(rho*c_p*dx)+(2*h*dt)/(rho*c_p*dy)) + T(i-1,j,k)*(2*a*dt)/(dy^2) + T(i,j+1,k)*(2*a*dt)/(dx^2);

% Upper Right Node (157)
        elseif j==n_j && i==1
           T(i,j,k+1)= T(i,j,k)*(1 - (2*h*dt)/(rho*c_p*dx)  - (2*a*dt)/(dx^2) - (2*a*dt)/(dy^2)) + T_inf*((2*h*dt)/(rho*c_p*dx)) + T(i+1,j,k)*(2*a*dt)/(dy^2) + T(i,j-1,k)*(2*a*dt)/(dx^2) + ((2*q_dot*dt)/(rho*c_p*dx));

% Lower Right Node (182)
        elseif j==n_j && i==n_i
           T(i,j,k+1)= T(i,j,k)*(1 - (2*h*dt)/(rho*c_p*dx) -(2*h*dt)/(rho*c_p*dy) - (2*a*dt)/(dx^2) - (2*a*dt)/(dy^2)) + T_inf*((2*h*dt)/(rho*c_p*dx)+(2*h*dt)/(rho*c_p*dy)) + T(i,j-1,k)*(2*a*dt)/(dx^2) + T(i-1,j,k)*(2*a*dt)/(dy^2) + ((2*q_dot*dt)/(rho*c_p*dx));

           end
       end
    end

    t = t+dt; %Sets the current time ahead 0.08 seconds until reaching 8 seconds
    k = k+1; %Tells the program to reiterate with the next time step

end

M = T(1:n_i,1:n_j,1:k); %Builds the 101 matrices for the temperature profile
    [T_max, Indx] = max(M(:)); %Returns the maximum value throughout all matrices
    [Indx1, Indx2, Indx3] = ind2sub(size(M), Indx); %Returns the location of T_max
        disp('Maximum Temperature:');
        disp(max(M(:)));
        disp('Location of Maximum Temperature:')
        disp([Indx1, Indx2, Indx3]);
FX = gradient(M);
    [FXmax, Indx] = max(FX(:));
    [Indx1, Indx2, Indx3] = ind2sub(size(FX), Indx);
        disp('Maximum Temperature Gradient:');
        disp(max(FX(:)));
        disp('Location of Maximum Temperature Gradient in x direction:')
        disp([Indx1, Indx2, Indx3]);

t = (k-1)*dt;

%% Plots

for i = 1:n_i
    y(i) = dy*(i-1) ;
end

for j = 1:n_j
    x(j) = dx*(j-1) ;
end

    set(gcf, 'units', 'normalized');
    set(gcf, 'Position', [0, 0.1, 1, 0.9]);


subplot(2,2,4)
pcolor(x,y,T(1:n_i,1:n_j,k))
    title('2D Temperature Profile')
    xlabel('Width (in)')
    ylabel('Height (in)')
    zlabel('Temperature (F)')
    set(gca,'XMinorTick','on')
    set(gca,'YMinorTick','on')
shading interp

subplot(2,2,1)
surf(x,y,T(1:n_i,1:n_j,k));
    title('3D Temperature Profile')
    xlabel('Width (in)')
    ylabel('Height (in)')
    zlabel('Temperature (F)')
    set(gca,'XMinorTick','on')
    set(gca,'YMinorTick','on')
    set(gca,'ZMinorTick','on')
    colorbar
shading interp

subplot(2,2,3)
[c,T] = contour(x,y,T(1:n_i,1:n_j,k));
    title('Temperature Contour')
    xlabel('Width (in)')
    ylabel('Height (in)')
    set(gca,'XMinorTick','on')
    set(gca,'YMinorTick','on')
    clabel(c,T);

subplot(2,2,2)