Crank Nicolson scheme using Finite Volume for 1D Transient Problem.

Dear Community,

I have attached discretized scheme (Problem_1_1 and Problem_1_2). According to complete MATLAB code which is working. I am getting suspicious about results when I compare them with books that are related to the Fully Implicit scheme. A few days earlier I posted my question. Any best idea which can do to satisfy?

Your help will be highly appreciated.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc
%% Sort out Inputs
Theta = 0.5; % Crank - Nicolson scheme
% Number of Control Volumes
N = 5;
% Domain length
L = 0.02; %[m]
% Grid Spacing (Wall thickness)
h = L/(N);  
% Diffusivity (Thermal conductivity)
k = 10;  %[W/m-K]
% Specific Heat Capacity
Row_c = 10e6;  %[J/m^(3)-K]
% The center of the first cell is at dx/2 & the last cell is at L-dx/2.
x = h/2 : h : L-(h/2);
% Diffusivity (Heat)
alpha = k/Row_c;
% Simulation Time
t_Final = 40;
t_steps = 8;
% Discrete Time Steps
dt = t_Final/(t_steps);
% Left Surface Temperature
dTdt = 0;    %[\circ c] 
% Right Surface Temperature
T_b = 0;    %[\circ c] 
% Parameteric Setup
lambda = (alpha.*dt)/(h^2)
ap_0 = Row_c * h/(dt);
dt_Check = Row_c *(h^2)/k;
% Initializing Variable
% C*T^(n+1) = BB*T^(n)
%
% At Node # 1 dpdx=0 is adjusted during simplification of scheme
%
% (1+lamda/2)T_P = (1/2)*lamda*T_E +(T_P)_0
%
% At Node # 5
%
% (1+2*lamda)T_P = lamda*T_W + lamda*T_b+(T_P)_0
%
% At Node Interior Nodes
%
% -(lamda/2)*T_W + (1+lamda)T_P -(lamda/2)*T_E = (lamda/2)*T_W_0 + (1-lamda)T_P_0 +(lamda/2)*T_E_0
%
% Unknowns at time level n
T_Old = zeros(N,1);
% Unknowns at time level n+1
T_New = zeros(N,1);
% Initial Value (Initial Condition)
T_Old(:) = 200; 
C = zeros(N,N);
D = zeros(N,1);
for t=1:dt:t_Final
for i = 1: N
   if  i > 1 && i < N
      C(i,i) = (1 + lambda);
      BB(i,i) = (1-lambda);
      C(i,i+1) = -0.5*lambda;
      BB(i,i+1) = 0.5*(lambda);
      C(i,i-1) = -0.5*lambda;
      BB(i,i-1) = 0.5*(lambda);
      D(i)=T_Old(i);
   elseif i == 1
    % For 1st boundary node
      C(i,i) = (1 + 0.5*lambda);
      C(i,i+1) = -0.5*lambda;
      D(i) = T_Old(i); 
     else
     % For Last boundary node
       C(i,i) = (1 + 2*lambda);
       C(i,i-1) = -lambda; 
       BB(i,i) = 0;
       D(i) = T_b*lambda + T_Old(i); 
   end
end
Temp = inv(C)*BB;
T_New = Temp*T_Old;
T_Old = T_New;
end 

T_New