I was wondering how to set up the backwards sweep (i=M:-1:1) to handle the values from the forward sweep. The matrix Tn1 updates to Tn, but then for the updated matrix after the backwards sweep the data isn't connected at all to the forward sweep.

function TwoDTriangularFinProject()

%% Geometry information

L= 0.10; %Length of Fin, meters

w= 0.20; %Width of Fin, meters

tt=0.01; %Full thickness of the Fin, meters

t= @(y)(tt/L)*y; %Thickness of Triangular Fin is a function of the y direction, meters

%% Property information

k= 200; %thermal conductivity, W/m-K

%% Boundary condition information

BCx0=1;%BC Flag at x=0,1 = constant temperature, 0= Specified Heat Transfer Coefficient(Convection), or Heat Flux

%Constant Convection Coeffient for Question 2 of Part 2 in Project, found analytically

BCy0=1;%BC Flag at y=0,1 = Constant temperature, 0= Specified Heat Transfer Coeffieient(Convection), or Heat Flux

%BCyL=0;%BC Flag at y=L,1 = Constant temperature, 0= Specified Heat Transfer Coeffieient(Convection), or Heat Flux

Tb= 120; %base temperature at y=0, deg C

%***In d vector for the tridiagonal matrix, any heat flux or convection terms

%present should be set to zero if there is an adiabatic BC!!***

%% Source term information

%h=@(x)7*x^(-1/2); %convection coefficient equation

h=[2213.5944,7.0000,4.9497,4.0415,3.5000,3.1305,2.8577,2.6458,2.4749,2.3333,2.2136,2.1106,2.0207,1.9415,1.8708,1.8074,1.7500,1.6977,1.6499,1.6059,1.5652,1.5275,1.4924,1.4596,1.4289,1.4000,1.3728,1.3472,1.3229,1.2999,1.2780,1.2572,1.2374,1.2185,1.2005,1.1832,1.1667,1.1508,1.1355,1.1209,1.1068]; %Calculated Convection Coefficient vector

Tinf= 30; %freestream temperature, deg C

%% Discretization information

M= 41; %number of nodes in the x direction

N= 31; %number of nodes in the y direction

dx=w/(M-1); %spacing between nodes, m

x=[0:dx:w];

dy=L/(N-1); %spacing between nodes, m

y=[0:dy:L];

%% Determine cross section and surface area for each CV

ds=zeros(length(y),1);

dAs=zeros(length(y),1);

Ac=zeros(length(y),1);

for j=2:N-1

t1=(t(y(j))+t(y(j+1)))/2; %average thickness of CV on side closer to y=0

t2=(t(y(j-1))+t(y(j)))/2; %avgerage thickness of CV on side closer to y=L

ds(j)= sqrt((L^.2)+((t(y(j-1)+t(y(j+1))))/2).^2)/L; %sqrt((t(y(j+1))-t(t(j-1))).^2+(2*dy).^2)/2; %using central difference

dAs(j)= (w+((t1+t2)/L))*ds(j)*dx; %surface area

Ac(j)=(t(y(j)))*dx; %Cross sectional area

end

%Surface Area and Area for the j = 1 location of CV (Different because it's area at the domains edge)

t1=t(y(1));

t2=(t(y(1))+t(y(2)))/2;

ds(1)=sqrt((L^2)+((t1+t2)).^2)/L;

dAs(1)=(w+((t1+t2)/L))*ds(1);

Ac(1)= tt*w; %(t(y(1)))*dx;

%Surface Area and Area for the j = N location of CV (Different because they area at the domains Tip)

t1=(t(y(end-1))+t(y(end)))/2;

t2=t(y(end));

ds(end)=sqrt((L.^2)+((t1+t2)).^2)/L;

dAs(end)=(w+((t1+t2)/L))*ds(end);

Ac(end)=(t(y(end)))*dx;

%% Start solution

Tn=Tinf.*ones(N,M); %initial guess (n) for domain temperature

Tn1=Tn; %initialize new iteration ("n+1") temperature

% Initialize tridiagonal vectors

a=zeros(N-1,1);

b=zeros(N,1);

c=zeros(N-1,1);

d=zeros(N,1);

% Loop through iterations, indexed by n

for n=1:100 %upper limit on iterations

% Sweep in positive x direction first, solving all j values

% simultaneously using tridiagonal solver

for i=1:M

%fill out first element of tridiagonal vectors

% a(1) = is absent because there is no j-1 term in the tridiagonal vector

b(1)= -((1/dy)*2*k*Ac(j+1)*k*Ac(1)/(k*Ac(j+1)+k*Ac(1))+(1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2*dy)))+(1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2*dy)))+((h(1))*dAs(j)));%Anorth+Awest+Aeast + Convection wrt x @base

c(1)= (1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2*dy)));%Aeast @base

d(1)=-((1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2))*dy)+(1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2))*dy)+((h(1))*dAs(j))); %Aeast+Awest + Convection wrt x @base

%fill out middle of tridiagonal vectors( **Boundary Nodes are already found from the boundary conditions, and the

%for loop stops at j = N-1 because there is a boundary at N, and since the area of the fin tip is basically nothing

%I'm going to assume that the temperature is equal to the free stream temperature directly at the fin tip**)

for j=2:N-1

a(j-1)= (1/dy)*2*k*Ac(j)*k*Ac(j-1)/(k*Ac(j)+k*Ac(j-1)); %Asouth, Harmonic mean of Areas

c(j)= (1/dy)*2*k*Ac(j+1)*k*Ac(j)/(k*Ac(j+1)+k*Ac(j)); %Anorth, harmonic mean of Areas

b(j)= -(a(j-1)+c(j)+(1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2)*dy))+(1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2)*dy))+((h(i))*dAs(j))); %Asouth+Anorth+Awest+Aeast + Convection wrt x

d(j)= -((1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2)*dy))+(1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2)*dy))+(h(i)*dAs(j)*Tinf)); %Aeast+Awest + Convection wrt x

end

%fill out last element of tridiagonal vectors (**Pretty Sure this is the

%Fin Tip Node FDE**)

%c(end) = is absent because there is no j+1 term in the tridiagonal vector

a(end)=(1/dx)*(k*(1/2*(t(end-1)+t(end))*dy));%Awest @Fin tip

b(end)= -((1/dy)*2*k*Ac(end)*k*Ac(end-1)/(k*Ac(end)+k*Ac(end-1))+(1/dx)*(k*(1/2*(t(end-1)+t(end))*dy))+(1/dx)*(k*(1/2*(t(end-1)+t(end))*dy))+((h(41))*dAs(end)*Tinf)); % Asouth+Awest+Aeast + Convection wrt x @Fin tip

d(end)= -((1/dx)*(k*(1/2*(t(end-1)+t(end))*dy))+(1/dx)*(k*(1/2*(t(end-1)+t(end))*dy))+((h(41))*dAs(end)*Tinf));% Aeast+Awest + Convection wrt x @fin tip

%solve tridiagonal matrix to update Tn+1 at current i location

Tn1(:,i)=tridiagonal(b,a,c,d);

end

%Initialize past solution at n with recently solved values at n+1 (** N+1 step is done

% after every i,j node is visited, it becomes the new Tn**)

Tn=Tn1;

% Sweep backward in negative x direction, solving all j values

% simultaneously using tridiagonal solver

for i=M:-1:1

%fill out first element of tridiagonal vectors

%c(end) = is absent because there is no j+1 term in the tridiagonal vector

a(end)=(1/dx)*(k*(1/2*(t(end-1)+t(end))*dy));%Awest @Fin tip

b(end)= -((1/dy)*2*k*Ac(end)*k*Ac(end-1)/(k*Ac(end)+k*Ac(end-1))+(1/dx)*(k*(1/2*(t(end-1)+t(end))*dy))+(1/dx)*(k*(1/2*(t(end-1)+t(end))*dy))+((h(41))*dAs(end)*Tinf)); % Asouth+Awest+Aeast + Convection wrt x @Fin tip

d(end)= -((1/dx)*(k*(1/2*(t(end-1)+t(end))*dy))+(1/dx)*(k*(1/2*(t(end-1)+t(end))*dy))+((h(41))*dAs(end)*Tinf));% Aeast+Awest + Convection wrt x @fin tip

%fill out middle of tridiagonal vectors

for j=2:N-1

a(j-1)= (1/dy)*2*k*Ac(j)*k*Ac(j-1)/(k*Ac(j)+k*Ac(j-1)); %Asouth, Harmonic mean of Areas

c(j)= (1/dy)*2*k*Ac(j+1)*k*Ac(j)/(k*Ac(j+1)+k*Ac(j)); %Anorth, harmonic mean of Areas

b(j)= -(a(j-1)+c(j)+(1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2)*dy))+(1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2)*dy))+((h(1))*dAs(j))); %Asouth+Anorth+Awest+Aeast + Convection wrt x

d(j)= -((1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2)*dy))+(1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2)*dy))+((h(1))*dAs(j)*Tinf)); %Aeast+Awest + Convection wrt x

end

%fill out last element of tridiagonal vectors

b(1)= -((1/dy)*2*k*Ac(j+1)*k*Ac(1)/(k*Ac(j+1)+k*Ac(1))+(1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2*dy)))+(1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2*dy)))+((h(1))*dAs(j)));%Anorth+Awest+Aeast + Convection wrt x @base

c(1)= (1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2*dy)));%Aeast @base

d(1)=-((1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2))*dy)+(1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2))*dy)+((h(1))*dAs(j))); %Aeast+Awest + Convection wrt x @base

%solve tridiagonal matrix to update Tn+1 at current i location

Tn1(:,i)=tridiagonal(b,a,c,d);

end

%Initialize past solution at n with recently solved values at n+1(** N+1 step is done

% after every i,j node is visited, it becomes the new Tn**)

Tn=Tn1;

%Check convergence of solution, etc

end

%% Appropriate BC's for vectors along y axis boundaries

if BCy0==1

b(1)=1;

c(1)=0;

d(1)=Tb;

else

b(1)= -((1/dy)*2*k*Ac(j+1)*k*Ac(1)/(k*Ac(j+1)+k*Ac(1))+(1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2*dy)))+(1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2*dy)))+((h(1))*dAs(j)));%Anorth+Awest+Aeast + Convection wrt x @base

c(1)= (1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2*dy)));%Aeast @base

d(1)=-((1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2))*dy)+(1/dx)*(k*(1/2*((t(j-1)+t(j))/2+(t(j)+t(j+1))/2))*dy)+((h(1))*dAs(j))); %Aeast+Awest + Convection wrt x @base

end

% if BCx0==1

% b(1)=0;

% c(1)=1;

% d(1)=Tinf;

% else

% %c(end) = is absent because there is no j+1 term in the tridiagonal vector

% a(end)=(1/dx)*(k*(1/2*(t(end-1/2)+t(end))*dy));%Awest @Fin tip

% b(end)= -((1/dy)*2*k*Ac(end)*k*Ac(end-1)/(k*Ac(end)+k*Ac(end-1))+(1/dx)*(k*(1/2*(t(end-1/2)+t(end))*dy))+(1/dx)*(k*(1/2*(t(end-1/2)+t(end))*dy))+(h(x(i))*dAs(end)*Tinf)); % Asouth+Awest+Aeast + Convection wrt x @Fin tip

% d(end)= -((1/dx)*(k*(1/2*(t(end-1)+t(end))*dy))+(1/dx)*(k*(1/2*(t(end)+t(end))*dy))+(h(x(i))*dAs(end)*Tinf));% Awest+Aeast + Convection wrt x @fin tip

end

% T = tridiagonal(b,a,c,d);

% figure;

% plot(h,'-b','LineWidth',2)

% xlabel('y [m]');

% ylabel('h(x) [w/m^2k]')

% ax=gca;

% ax.FontSize=14;

%Computing the heat flux

%Displaying the results

function x = tridiagonal(diag,lower,upper,rhs)

% tridiagonal: solves a tridiagonal matrix.

%

% Synopsis: x = tridiagonal(diag,lower,upper,rhs)

%

% Input: diag = diagonal vector of length N

% lower = lower diagonal vector of length N-1

% upper = upper diagonal vector of length N-1

% rhs = right hand side vector of length N

%

% Output: x = solution vector of length N

%

% --- determine the number of equations

N = length(diag);

% --- eliminate the lower diagonal

for i = 2:N

factor = lower(i-1)/diag(i-1);

diag(i) = diag(i)-upper(i-1)*factor;

rhs(i) = rhs(i)-rhs(i-1)*factor;

end

% --- back substitute

x(N) = rhs(N)/diag(N);

for i = N-1:-1:1

x(i) = ( rhs(i)-upper(i)*x(i+1) ) / diag(i);

end

end