While i agree with Wolfgang that its best to choose Crank-Nicolson time stepping, i think it is incorrect to assert explicit time stepping results in 'incredibly small' time steps or even to suggest your script is not functioning because of explicit time stepping. You can make it work just fine, but you need to get your discretization correct and to implement the boundary conditions properly.
Problems:
First, your discretization of your spatial second derivative is incorrect. I like to do this using the Finite Volume method because you can not go wrong. You define an average temperature over a control volume $[x-\frac{1}{2}\Delta x,x+\frac{1}{2}\Delta x]$ where $x$ is the node at the center of the volume element with size $\Delta x$:
$$\bar{u} = \frac{1}{\Delta x}\int_{x-\frac{1}{2}\Delta x}^{x+\frac{1}{2}\Delta x} u dx$$
Now you integrate the heat equation:
$$\partial_{t}u=a\partial_{x}^{2}u$$
over the control volume:
$$\frac{1}{\Delta x}\int_{x-\frac{1}{2}\Delta x}^{x+\frac{1}{2}\Delta x} \partial_{t}u dx=a\frac{1}{\Delta x}\int_{x-\frac{1}{2}\Delta x}^{x+\frac{1}{2}\Delta x} \partial_{x}^{2}u dx$$
$$\partial_{t}\bar{u} =\frac{a}{\Delta x}\left.\partial_{x}u \right|_{x-\frac{1}{2}\Delta x}^{x+\frac{1}{2}\Delta x}$$
Edit: The result here is that the temperature in the control volume is spatially uniform, i.e. the average temperature $\bar{u}$ at node $x$. For $\Delta x$ small enough, this become increasingly more accurate. The change in time of $\bar{u}$ then becomes dependent of the incoming and outgoing flux of heat at the boundaries located at $x-\frac{1}{2}\Delta x$ and $x+\frac{1}{2}\Delta x$. The choice of these coordinates make this a staggered grid where the nodes do not align with the boundaries.
using the second-order accurate discretization of the derivate, we get:
$$\partial_{t}\bar{u}
= \frac{a}{\Delta x}\left[\frac{\bar{u}\left(x+\Delta x\right)-\bar{u}\left(x\right)}{\Delta x}-\frac{\bar{u}\left(x\right)-\bar{u}\left(x-\Delta x\right)}{\Delta x}\right]
= \frac{a}{\Delta x^{2}}\left[\bar{u}\left(x+\Delta x\right)-2\bar{u}\left(x\right)+\bar{u}\left(x-\Delta x\right)\right]
$$
which is different from your discretization. Edit: since the averaged temperature at node $x$ is used exclusively in these equations i will drop the 'bar' from $\bar{u}$
Edit: For abritrary coordinate $\tilde{x}$, second-order accuracy is shown by Taylor expansion:
$$u\left(\tilde{x}-\frac{\Delta x}{2}\right)=u\left(\tilde{x}\right)-\frac{\Delta x}{2}\left.\frac{du}{dx}\right|_{\tilde{x}}+\frac{1}{2}\left(\frac{\Delta x}{2}\right)^{2}\left.\frac{d^{2}u}{dx^{2}}\right|_{\tilde{x}}+O\left(\Delta x\right)^{3}$$
$$u\left(\tilde{x}+\frac{\Delta x}{2}\right)=u\left(\tilde{x}\right)+\frac{\Delta x}{2}\left.\frac{du}{dx}\right|_{\tilde{x}}+\frac{1}{2}\left(\frac{\Delta x}{2}\right)^{2}\left.\frac{d^{2}u}{dx^{2}}\right|_{\tilde{x}}+O\left(\Delta x\right)^{3}$$
Subtracting to cancel the second derivative gives:
$$u\left(\tilde{x}+\frac{\Delta x}{2}\right)-u\left(\tilde{x}-\frac{\Delta x}{2}\right)=\Delta x\left.\frac{du}{dx}\right|_{\tilde{x}}+O\left(\Delta x\right)^{3}$$
rewriting to an expression for the first derivative:
$$\left.\frac{du}{dx}\right|_{\tilde{x}}=\frac{u\left(\tilde{x}+\frac{\Delta x}{2}\right)-u\left(\tilde{x}-\frac{\Delta x}{2}\right)}{\Delta x}+O\left(\Delta x\right)^{2}$$
substituting in the averaged heat equation and filling in the $\tilde{x}=x-\frac{1}{2}\Delta x$ and $\tilde{x}=x+\frac{1}{2}\Delta x$ respectively gives the discretized equation found before.
Secondly, you need to define boundary conditions. Since you didn't specify i am going to assume you assume zero-gradient boundary conditions. Since $x$ is at the center of the node, you can view it as a staggered grid. This means that the boundaries of the domain are in between nodes at $x_{b0}=x_1-\frac{1}{2}\Delta x$ and $x_{bf}=x_N+\frac{1}{2}\Delta x$, where $x_{1}$ and $x_{N}$ are the first and $N$th node respectively. We can set the gradients at $x_{b0}$ and $x_{bf}$ by introducing 'ghost' nodes outside the domain at $x_0$ and $x_{N+1}$ and then the gradients are defined as:
$$\partial_{x}u\left(x_{b0}\right)=\frac{u\left(x_{1}\right)-u\left(x_{0}\right)}{\Delta x}=0\quad\partial_{x}u\left(x_{bf}\right)=\frac{u\left(x_{N+1}\right)-u\left(x_{N}\right)}{\Delta x}=0$$
So simply setting $u\left(x_{0}\right)=u\left(x_{1}\right)$ and $u\left(x_{N+1}\right)=u\left(x_{N}\right)$ will ensure zero-gradients at the boundaries.
Now we are in a position to solve the equations. Edit: For each node we solve the discretized heat equation and use the information at the boundaries for the first and last node. For simplicity each node is identified by a subscript $i$, i.e. $u_i=u\left(x_i\right)$:
$$\partial_{t}\bar{u}_{i}=\frac{a}{\Delta x^{2}}\left[\bar{u}_{i+1}-2\bar{u}_{i}+\bar{u}_{i-1}\right]$$
written out explicitly:
$$\partial_{t}\bar{u}_{1} = \frac{a}{\Delta x^{2}}\left[\bar{u}_{2}-2\bar{u}_{1}+\bar{u}_{0}\right]$$
$$\partial_{t}\bar{u}_{2} = \frac{a}{\Delta x^{2}}\left[\bar{u}_{3}-2\bar{u}_{2}+\bar{u}_{1}\right]$$
$$\vdots $$
$$\partial_{t}\bar{u}_{N-1} = \frac{a}{\Delta x^{2}}\left[\bar{u}_{N}-2\bar{u}_{N-1}+\bar{u}_{N}\right]$$
$$\partial_{t}\bar{u}_{N} = \frac{a}{\Delta x^{2}}\left[\bar{u}_{N+1}-2\bar{u}_{N}+\bar{u}_{N-1}\right]$$
Edit: Here $u_0$ and $u_{N+1}$ are the values of the temperature at the ghost nodes previously determined from the boundary conditions. We substitute those in to get:
$$\partial_{t}\bar{u}_{1} = \frac{a}{\Delta x^{2}}\left[\bar{u}_{2}-\bar{u}_{1}\right]$$
$$\partial_{t}\bar{u}_{N} = \frac{a}{\Delta x^{2}}\left[-\bar{u}_{N}+\bar{u}_{N-1}\right]$$
Next, we define a system of equations with vector:
$$\vec{u}=\begin{bmatrix}u\left(x_{1}\right) & u\left(x_{2}\right) & \cdots & u\left(x_{N-1}\right) & u\left(x_{N}\right)\end{bmatrix}^{T}$$
which allows us to write these equations as a system of equations:
$$\partial_{t}\begin{bmatrix}\bar{u}_{1}\\
\bar{u}_{2}\\
\vdots\\
\bar{u}_{N-1}\\
\bar{u}_{N}
\end{bmatrix}=\frac{a}{\Delta x^{2}}\begin{bmatrix}-1 & 1\\
1 & -2 & 1\\
& \ddots & \ddots & \ddots\\
& & 1 & -2 & 1\\
& & & 1 & -1
\end{bmatrix}\begin{bmatrix}\bar{u}_{1}\\
\bar{u}_{2}\\
\vdots\\
\bar{u}_{N-1}\\
\bar{u}_{N}
\end{bmatrix}$$
so we get the linear system:
$$\partial_t\vec{u}=\boldsymbol{A}\cdot\vec{u}$$
where:
$$A=\frac{a}{\Delta x^{2}}\left[\begin{array}{ccccc}
-1 & 1\\
1 & -2 & 1\\
& \ddots & \ddots & \ddots\\
& & 1 & -2 & 1\\
& & & 1 & -1
\end{array}\right]$$
All rows are described by our discretized equation, except in the first and last column where the values for the 'ghost' nodes outside of the domain are used to account for the boundary conditions.
Let's not worry about discretization of the time derivative, we will simply let matlab integrate this system using an appropriate 'ode' function. This has as advantage that matlab choses an proper timestep for discretization which makes sure that the integration is stable.
Below you will find an implementation of the above linear system and the result.
Code and result:
function main()
clc, clear all, close all;
a = 1; %% thermal diffusivity
pt=100; %%% number of points
xsta=-5; %%% x start
xend=5; %%% x end
x = linspace(xsta,xend,pt); %%% creating x space -5 ~ 5
% create discretization matrix
A = zeros(pt,pt);
A(1,1) = -1;
A(1,2) = 1;
for i = 2:pt-1
A(i,i-1) = 1;
A(i,i) = -2;
A(i,i+1) = 1;
end
A(pt,pt-1) = 1;
A(pt,pt) = -1;
% initial condition
u0 = zeros(pt,1);
u0(abs(x)<1) = 1;
function dudt = pde(t, u, Fo_x)
dudt = Fo_x*A*u;
end
nt = 5; % increase nt if you want more time points (not needed for accuracy)
Fo = 1; % Fo is a dimensionless diffusion coefficient, a*t/l^2
Fo_x = Fo*(pt-1); % Fo_x is diffusion coefficient based on grid
tspan = linspace(0,1,nt);
[dimt,u] = ode45(@pde, tspan, u0, [], Fo_x);
figure, plot(x,u,'o')
legendCell = cellstr(num2str((dimt*(xend-xsta)^2/a), 't=%-d s'));
legend(legendCell)
xlabel('spatial coordinate, x')
ylabel('temperature, u')
end
