Help solving a PDE using the Thomas algorithm (tridiagonal matrix).

So I'm writing a script to solve for the following 1-D unsteady heat equation with the given boundary and initial conditions:

For the numerical simulation, it's required to use the following combined method:

To solve the solution above, it's necessary to use the Thomas algorithm.

Below shows the variable values during different "cases"

Here's the code I have so far. It gives me a matrix for the temperature T with only the initial condition and boundaries filled, the rest are zero. I need to solve for the remaining values.

Please help!

% Inputs
tend = 30; %time end 
L = 10;
alpha = 2;

% Cases:
R = [1, 1, 1, 1, 1/2, 1/(sqrt(20))];
theta = [1, 3/4, 1/2, 5/12 ,1/3,(1/2)-(1/(12*R(6)))];
nx = 101;

for k = 1 % for 6 cases
  r = R(k);
  th = theta(k);
  dx = L/(nx-1);
  dt = r*dx^2/alpha;
  x = (0:dx:L);
  t = (0:dt:tend);
  nt = length(t);

% Initialization
  T = zeros(nx,nt);

% inital conditions
  T(:, 1) = 100 * sin(pi*x/L);

% boundary conditions
  T(1, :) = 0;
  T(nx,:) = 0;

   ...
   . 
   . 
   end