Finite differences for the 1-d heat equation

c:1.2;

/* subdivisions */
L:.6;
n:18;
dx:L/n;

dt:0.0005;

h:dt/(dx)^2;

A:genmatrix(lambda([i,j],
if (i=1 or i=n+1) then
    (
    if (j=i)
    then 1
    else 0
    )
else
    (
    if (j=i)
    then 1-2*c^2*h
    else
        (
        if (j=i+1 or j=i-1) then c^2*h
            else 0
        )
    )
),n+1,n+1)$

 /* frequency multiplier */
k:2;

/* pure frequency exact solution with BC u(0,t)=0, u(L,t)=1 */
u(x,t):=sin(%pi*k*x/L)*exp(-(c*%pi*k/L)^2*t)+x/L;
/* plot2d([parametric,x,u(x,0),[x,0,L]])$ */

xs:create_list(i*dx,i,0,n)$
us:u(xs,0)$float(%)$
/* plot2d([[discrete,xs,us]])$ */

plt:[[discrete,xs,us]]$
/* plot2d(plt)$ */

/* T = time step to plot */
T:2;
for i:1 thru T do
    (
    us:A.us,
    if (i=T) then plt:append(plt,[[discrete,xs,create_list(us[i][1],i,1,n+1)]])
    )$

/* exact solutions */
plt:append(plt,[[parametric,x,u(x,0),[x,0,L]]])$
plt:append(plt,[[parametric,x,u(x,T*dt),[x,0,L]]])$


plot2d(plt)$