How can I solve my 1D finite volume heat transfer problem using the Crank-Nicolson method?

Question

I was hoping to get some help solving the heat transfer equation using the Crank-Nicolson method.

I have my heat transfer equation in linear form:

Equation 1

\frac{\partial T}{dt}=\omega \left ( A\cdot T+b \right )

I am able to solve this directly using ode15s in matlab. Howerver for my application I need to be able to do this within a for loop with a prescribed time increment. I am trying to insert this into the following equation in order to solve it using the Crank Nicholson implicit method:

Equation 2

q^{n+1}=q^{n}+\Delta t\left [ \phi f\left ( t_{n+1},u^{n+1} \right )+\left ( 1-\phi \right )f\left ( t_{n},u^{n} \right ) \right ] where \phi =\frac{1}{2}

However the results produced appear to oscillate around the true temperature values. Hope somebody can help! Thank you in advance.

%% Single material with internal heat generation
clc
clear

% 1. SET UP PROBLEM DIMENSIONS
L = 1; % Length in m
N = 100; % number of nodes
h = L/N; % node length

% 2. SET UP TEMPERATURES AT x=0 and x=L (IF USING THERMAL BCs)
TA = 1;
TB = 1;

% 3. SET UP MATERIAL ZONES
nZones = 2;
zoneBoundaries = [1,50;51,100];
zoneMap = zeros(N,1);
for i = 1:nZones
zoneMap(zoneBoundaries(i,1):zoneBoundaries(i,2)) = i;
end

% 4. SET UP MATERIAL PROPERTIES
kWood = 10; % Thermal conductivity of Wood / W m^-1 K^-1
rhoWood = 1000; % Density of Wood [kg/m3]
CpWood = 100; % Heat capacity of Wood [J/kgK]

kSteel = 45; % Thermal conductivity of steel / W m^-1 K^-1
rhoSteel = 7850; % Density of steel [kg/m3]
CpSteel = 420; % Heat capacity of steel [J/kgK]

% 5. SET UP THE MATRICES FOR HEAT EQUATION
A = zeros(N, N);
A(1, 1) = -3;
A(1, 2) = 1;
A(N, N) = -3;
A(N, N-1) = 1;
for i = 2:N-1
    A(i, i - 1) = 1;
    A(i, i) = -2;
    A(i, i + 1) = 1;
end

% Adjust A for material boundaries
kappa = kWood/kSteel;
gamma = (kappa-1) / (kappa+1);

A(50,:) = 0;
A(50,50-1) = 1;
A(50,50) = -1*(2-gamma);
A(50,50+1) = 2/(kappa+1);

A(51,:) = 0;
A(51,51-1) = 2*kappa/(kappa+1);
A(51,51) = -1*(2+gamma);
A(51,51+1) = 1;

b = zeros(N, 1);
b(1) = 2*TA;
b(end) = 2*TB;

% Set up omega matrix
k = zeros(N, 1);
k(zoneMap == 1) = kWood;
k(zoneMap == 2) = kSteel;
rho = zeros(N, 1);
rho(zoneMap == 1) = rhoWood;
rho(zoneMap == 2) = rhoSteel;
Cp = zeros(N, 1);
Cp(zoneMap == 1)  = CpWood;
Cp(zoneMap == 2)  = CpSteel;
alpha =  k./(rho.*Cp);

omega1 = zeros(N, N);
for i = 1:N
    omega1(i,i) = alpha(i);
end
omega = (1/h^2) .* omega1;

S = zeros(N, 1); % Rate of heat generation per kg J s^-1 kg^-1

% Differential Equation Function
diff_eq = @(t, q) omega * (A * q + b) + S;
%% SOLVE diff_eq ode15s (working)
t = logspace(-4, 3, 12);
T0 = zeros(N, 1);

sol = ode15s(diff_eq, [0 max(t)], T0);
i = linspace(0, N-1, N);
x = L / N * (i + 1 / 2);
yPlot = deval(sol, t);

colors = turbo(12);

figure
for i = 1:length(t)
    plot(x, yPlot(:,i), '-', 'DisplayName', ['t = ' num2str(t(i)) ' s'], 'Color', colors(i,:));
    hold on
end
legend('box', 'off')
xlabel('Spatial coordinate, x');
ylabel('Temperature coordinate, T');

%% SOLVE diff_eq Crank Nicolson (not working)
tStart = 0;
tEnd = 1000; % Set your desired end time
dt = 1; % Set your desired time step

% Discretize time
time = tStart:dt:tEnd;
numTimeSteps = length(time);

% Initial temperature
T0 = zeros(N, 1);

% Initialize solution matrix q at each time step
qMatrix = zeros(N, numTimeSteps);
qMatrix(:,1) = T0;

% Time-stepping loop
for i = 2:numTimeSteps
qMatrix(:,i) =  (1-0.5*dt*omega*A) \ ( (1 + 0.5*dt*omega*A)*qMatrix(:,i-1) + dt*omega*b + dt*S );
end

ii = linspace(0, N-1, N);
x = L / N * (ii + 1 / 2);

colors = turbo(length(1:size(qMatrix,2)));
figure
for i = 1:size(qMatrix,2)
    plot(x, qMatrix(:,i), '-','Color', colors(i,:));
    hold on
end
xlabel('Spatial coordinate, x');
ylabel('Temperature coordinate, T');

The equation produces sensible results when solved using matlabs built in ode solvers, however my implementation of the Crank-Nicolson method doesn't seem to work. The answer looks as though it is oscilating around the correct answer from one timestep to the next.