Here is the Matlab code for the heat and wave equations:

%This function loads EM2 into the sparse matrix EM2
function EM2=loadEM2

x_string='EM2_x.dat';
i_string='EM2_i.dat';
j_string='EM2_j.dat';

load 'EM2_x.dat';
load 'EM2_i.dat';
load 'EM2_j.dat';

EM2=sparse(EM2_i,EM2_j,EM2_x);

%This function loads IP2 into the sparse matrix IP2...

function IP2=loadIP2

x_string='IP2_x.dat';
i_string='IP2_i.dat';
j_string='IP2_j.dat';

load 'IP2_x.dat';
load 'IP2_i.dat';
load 'IP2_j.dat';

IP2=sparse(IP2_i,IP2_j,IP2_x);

%This function generates time series of a solution to the heat equation %k is the number of iterations you want to use (t=hk, where h is the %time interval, and k is the number of intervals). Make it bigger for %better but slower results. EM is the energy matrix, IP is the inner %product matrix, and c0 is the inital values. timestep is the %interval of time between steps of the time series, and numsteps is the %total number of steps. The results are stored in a matrix whose rows %are the splines representations of the function at a particular time.

function A=generate_heat_time_series(EM,IP,c0,k,timestep,numsteps)

A=[c0']

for i=1:numsteps
outC=heat_eqn_t(EM,IP,c0,timestep*i,k);
A=[A ; outC'];
end

%This function is just like generate_heat_time_series, except that it generates %the wave equation time series. All the parameters work in the same way, except that %you need to pass c0 and d0, with c0 being the initial values, and d0 being the %initial derivatives. For documentation on k, timestep and numsteps, see %generate_heat_time_series.m

function A=generate_wave_time_series(EM,IP,c0,d0,k,timestep,numsteps)

A=[c0']

for i=1:numsteps
outC=wave_eqn_t(EM,IP,c0,d0,timestep*i,k);
A=[A ; outC'];
end

%This function is just like generate_wave_time_series.m, except that %it calls wave_eqn_t2 instead of wave_eqn_t, which means that instead of %using the iterative type solutions, it uses the matrix exponential solutions %This has the advantage of being more accurate (in particular, you don't have to %worry about balancing space and time intervals) with the added cost of being %somewhat slower. Notice that the parameter k here is meaningless, since k is %not used when using the matrix exponential command...

function A=generate_wave_time_series2(EM,IP,c0,d0,k,timestep,numsteps)

A=[c0']

for i=1:numsteps
outC=wave_eqn_t2(EM,IP,c0,d0,timestep*i,k);
A=[A ; outC'];
end

%This function returns the HS norm of a matrix

function n=HSnorm(M)

thesize=size(M);
dimen=thesize(1);
n=0;

for i=1:dimen
for j=1:dimen
n=n+(M(i,j))^2;
end
end

n=sqrt(n);

%This function returns the spectral radius of a matrix M by taking %HSnorm(M^n)^1/n. The actual spectral radius is the limit by taking n=infinity. %However, it seems to converge from above, so you will be able to bound it by using %n=20 or something. Note that one must be careful using this function, because %depending on the values of n and M, you might get infinities and the such...

function rad=spec_radius(M,n)

rad=HSnorm(M^n)^(1/n);

%This function solves the heat eqn at time t, with t=hk. Here, h is the time interval, %and k is the number of time intervals. c0 is the intial values of the function...

function c=heat_eqn(EM,IP,c0,h,k)

thesize=size(EM);
dimen=thesize(1);

M=eye(dimen)-h*inv(IP)*EM;

c=(M^k)*c0;

%This function is just like heat_eqn.m, except that it takes t and k as parameters, %rather than h and k, thus allowing you to specify a time directly.

function c=heat_eqn_t(EM,IP,c0,t,k)

h=t/k;

thesize=size(EM);
dimen=thesize(1);

M=eye(dimen)-h*inv(IP)*EM;

c=(M^k)*c0;

%This function solves the wave eqn at time t with accuracy determined by h=t/k. %It takes c0 as the initial values, and d0 as the intial velocities.

function c=wave_eqn_t(EM,IP,c0,d0,t,k)

h=t/k;

thesize=size(EM);
dimen=thesize(1);

M=eye(dimen*2)+h*[zeros(dimen) eye(dimen);-1*inv(IP)*EM zeros(dimen)];

c_d=(M^k)*[c0;d0];

c=c_d(1:dimen);

%This function solves the wave eqn at a particular time t, just like %wave_eqn_t, except that instead of using k and thus h, it uses matrix %exponential. The advantage is that you don't have to worry about coordinating %the size of the time interval with the space size (size of the matrix), because %the matrix exponential will automatically go out to the required level of accuracy. %However, for the same reason, it is slower.

function c=wave_eqn_t2(EM,IP,c0,d0,t,k)

h=t/k;

thesize=size(EM);
dimen=thesize(1);

M=t*[zeros(dimen) eye(dimen);-1*inv(IP)*EM zeros(dimen)];

c_d=expm(full(M))*[c0;d0];

c=c_d(1:dimen);