% expm9.m - for Table 6
% This computes \norm T_t1-1\norm_\infty
% by discretizing the interval (0,a)



clear all




% Set up basic data

%h=number of divisions per unit length 
h=20

%a=length of interval
a=10

%b=diffusion paramter
b=20

% number of points altogether
N = h*a-1
 
 
 
 
 
% set up tridiagonal matrix

  u=ones(N,1);
  v=ones(N-1,1);
  A=-2*(h^2)*b^(-1)*diag(u)+(h^2)*b^(-1)*diag(v,1)+(h^2)*b^(-1)*diag(v,-1);
  B=(-h/2)*diag(v,1)+(h/2)*diag(v,-1);
  C=A+B;
  
% compute the quantity of interest
  
  format long
   rho=zeros(11,1);
   for tr=1:11     
      t=(tr-1)/20;
      rho(tr,1)=norm(expm(t*C)*ones(N,1)-ones(N,1),inf);
   end;
rho
   
hold on;
plot((0:10)/20,rho);