How does Matlab handle the edge of an array?

I have some code to preform the Crank-Nicolson Method on on arrays in both 1D and 2D, and while I was testing it for the case for the heat equation (thermal diffusion) I noticed some behaviour that didn't match the physics. This got me thinking about how Matlab handles the edges of arrays.

With thermal diffusion, if you start with one half of your array/grid at one uniform temperature and the other and a lower uniform temperature you expect the whole grid to reach some uniform temperature between the two. However I have noticed my code always reaches a global temperature of 0 for long enough times. In fact if I start with a grid all at a uniform temperature I still reach a global zero if left long enough. Almost as if heat is leaking off the sides of my array!

I have attached my code below in case my problem just down to how I have coded this. However could someone please explain how Matlab deals with the edges of 1 and 2D arrays. Am I correct in thinking Matlab will consider anything off my array to be a zero by default and is there a way of 'isolating' and array for cases like this (i.e. the heat equation).

%%Heat Equation Test
clear;                      
len = 2^6; lim = 0.1;       % Meters
dx = 2*lim/len; dt = 1;     % Meters; Seconds
steps = 100;                % Total Time Steps
x = linspace(-lim,lim,len);

% -------------------------------------------------------------------------     
u = zeros(len,1);
for i = 1:len
    if x(i) > (-lim/2) && x(i) < (lim/2)
        u(i) = 1; % Rectangular function (width a, centred about origin)
    end
end
% -------------------------------------------------------------------------
% Tridiagonal Matrix Elements:
a = 1.27*10^(-4);       % Gold's thermal diffusivity 
r = (a*dt)/(dx^2);      % Constant arising from Crank Nicolson Method

Aon = 2*(1+r); Aoff = -r;
Bon = 2*(1-r); Boff = r;

% Create Sparse Tridiagonal Matrices 
D = sparse(1:len,1:len,Aon*ones(1,len),len,len);
E = sparse(2:len,1:len-1,Aoff*ones(1,len-1),len,len);
A = E+D+E';

D = sparse(1:len,1:len,Bon*ones(1,len),len,len);
E = sparse(2:len,1:len-1,Boff*ones(1,len-1),len,len);
B = E+D+E';

U = B*u;
Unext = A\U;
for i = 1:steps  
    U = B*Unext;
    Unext = A\U;
end

Here is an example of the output of this code:

Physically you would obviously expect to get a flat line at 0.5 for lone enough times. If anyone can spot why my code doesn't give me this physical result I would greatly appreciate it!