SUPPLEMENTARY MATERIALS -- 14 MATLAB/CHEBFUN CODES
            Jared L. Aurentz and Lloyd N. Trefethen
         "Block operators and spectral discretizations"
                 SIAM Review, Education section
                         August 2016

% abbreviations.m - CONVENIENT STARTUP COMMANDS
MS = 'markersize'; FS = 'fontsize'; IN = 'interpreter';
LT = 'latex'; LW = 'linewidth';
fs = 13;
subplot(2,1,1)

% plottitle.m - PUTS TITLES ON PLOTS WITH "TEXT"
function plottitle(n,fs,pos)
aa = axis;
a = aa(1); b = aa(2); c = aa(3); d = aa(4);
if nargin==2  % title at right
  xpos = a+.97*diff([a b]); ypos = c+.88*diff([c d]); 
  text(xpos,ypos,['Example ' int2str(n)],...
    'fontsize',fs,'fontname','helvetica','horizontalalignment','right')
else          % title at left
  xpos = a+.03*diff([a b]); ypos = c+.88*diff([c d]); 
  text(xpos,ypos,['Example ' int2str(n)],...
    'fontsize',fs,'fontname','helvetica','horizontalalignment','left')
end
  
% ex01.m - EXAMPLE 1
abbreviations
f = @(x) exp(x); alpha = 2; beta = 0; n = 18;
L = 0.2*diffmat([n n+2],2) + diffmat([n n+2],1) + diffmat([n n+2],0);
vT = diffrow(n+2,0,-1); wT = diffrow(n+2,1,1);
A = [L; vT; wT];
rhs = [gridsample(f,n); alpha; beta];
u = A\rhs;
plot(chebfun(u),'.-',LW,.7,MS,16)
xlabel('$x$',FS,16,IN,LT), ylabel('$u$',FS,16,IN,LT), set(gca,FS,11)
plottitle(1,fs), print -depsc ex01

% ex02.m - EXAMPLE 2
abbreviations
n = 85; a = -20; b = 10; X = [a,b];
x = @(x) x;
L = diffmat([n n+2],2,X) - ...
   spdiags(gridsample(x,n,X),0,n,n)*diffmat([n n+2],0,X);
vT = introw(n+2,X);
wT = diffrow(n+2,0,b,X) - diffrow(n+2,0,a,X);
A = [L; vT; wT];
rhs = [zeros(n,1); 1; 0];
u = A\rhs;
plot(chebfun(u,X),'.-',LW,.7,MS,14)
xlabel('$x$',FS,16,IN,LT), ylabel('$u$',FS,16,IN,LT)
set(gca,FS,11), ylim([-.6 .6]), set(gca,'ytick',-.6:.3:.6)
plottitle(2,fs), print -depsc ex02

% ex03.m - EXAMPLE 3
abbreviations
n = 80; d = 10; X = [-d,d]; x = chebfun('x',X);
L = -diffmat([n n+2],2,X) + ...
     spdiags(gridsample(x.^2,n,X),0,n,n)*diffmat([n n+2],0,X);
vT = diffrow(n+2,0,-d,X); wT = diffrow(n+2,0,d,X);
A = [L; vT; wT];
I = diffmat([n n+2],0,X);
B = [I; zeros(2,n+2)];
[V,D] = eig(A,B);
[vals,ii] = sort(diag(D));
V = V(:,ii); D = D(ii,ii);
plot(x.^2,'-k',LW,3), hold on
for j = 1:6
  v = V(:,j);
  if v(end-30)<0, v = -v; end
  u = .6*v + D(j,j);
  plot(chebfun(u,X),'.-',LW,.7,MS,7)
  disp(D(j,j))
end
ylim([0 13.7]), set(gca,FS,11), set(gca,'ytick',0:4:12), hold off
xlabel('$x$',FS,16,IN,LT), ylabel('$u$',FS,16,IN,LT)
plottitle(3,fs), print -depsc ex03

% ex04.m - EXAMPLE 4
abbreviations
n = 50; d = 10; X = [-d,d]; x = chebfun('x',X);
A = -diffmat([n n],2,X,'trig') + diag(gridsample(x.^2,n,X,'trig'));
[V,D] = eig(A);
[vals,ii] = sort(diag(D));
V = V(:,ii); D = D(ii,ii);
plot(x.^2,'-k',LW,3), hold on
for j = 1:6
  v = V(:,j);
  if v(end-15)<0, v = -v; end
  disp(D(j,j))
  u = 1.6*v + D(j,j);
  plot(chebfun(u,X,'trig'),'.-',LW,.7,MS,7)
end
ylim([0 13.7]), set(gca,FS,11), set(gca,'ytick',0:4:12), hold off
xlabel('$x$',FS,16,IN,LT), ylabel('$u$',FS,16,IN,LT)
plottitle(4,fs), print -depsc ex04

% ex05.m - EXAMPLE 5
abbreviations
n = 90;
D = diffmat([n n+1]);
vT = diffrow(n+1,0,-1); wT = diffrow(n+1,0, 1); zT = 0*vT;
I = diffmat([n n+1],0);
x = chebfun('x');
XI = diag(gridsample(exp(3*x),n))*I;
A = [D -10*I; 30*XI D; vT zT; wT zT];
rhs = [zeros(2*n+1,1); 1];
uv = A\rhs;
u = uv(1:n+1); v = uv(n+2:2*n+2);
plot(chebfun(u),'.-k',LW,.7,MS,12), hold on
plot(chebfun(v),'.-',LW,.7,MS,12,'color',[0 .7 .2])
xlabel('$x$',FS,16,IN,LT), ylabel('$u$ and $v$',FS,16,IN,LT)
set(gca,FS,12), hold off
plottitle(5,fs,'left'), print -depsc ex05

% ex06.m - EXAMPLE 6  [there is no code for this]

% ex07.m - EXAMPLE 7
abbreviations
f = @(x) exp(x);
alpha = 2; beta = 0; n = 18;
x = chebfun('x');
L = 0.2*intmat(n,0) + intmat(n,1) + intmat(n,2);
vT = zeros(1,n); wT = introw(n);
A = [ [L; vT; wT] [gridsample(x+2,n); 0; 1] [ones(n,1); 1; 0] ];
rhs = [gridsample(f,n); alpha; beta];
vcd = A\rhs;
v = vcd(1:n); c = vcd(n+1); d = vcd(n+2);
vc = chebfun(v);
uc = cumsum(cumsum(vc)) + c*(1+x) + d;
plot(uc,'.-',LW,.7,MS,16)
xlabel('$x$',FS,16,IN,LT), ylabel('$u$',FS,16,IN,LT)
set(gca,FS,11,'xtick',-1:.5:1)
plottitle(7,fs), print -depsc ex07

% ex08.m - EXAMPLE 8
abbreviations
f = @(x) 1 + 0*x; g = @(x) 1 + 0*x;
w = @(x) exp(-3*x.^2); wf = @(x) f(x).*w(x);
c = 0; n = 25;
gT = introw(n);
A = [diag(gridsample(w,n)) gridsample(g,n) ; gT 0];
rhs = [gridsample(w,n); c]; ulam = A\rhs;
u = ulam(1:n);
plot(chebfun(u),'.-',LW,.7,MS,16), hold on
plot(chebfun(f),'--r',LW,1)
xlabel('$x$',FS,16,IN,LT), ylabel('$u$',FS,16,IN,LT)
set(gca,FS,11,'ytick',-4:2:2), ylim([-4 2.5]), hold off
plottitle(8,fs), print -depsc ex08

% ex09.m - EXAMPLE 9
abbreviations
alpha = 1; beta = 4/3; n = 20;
x = chebpts(n+2); f = exp(2*chebpts(n));
D2 = diffmat([n n+2],2); I = diffmat([n n+2],0);
N = @(u) D2*u + (I*u).^2 - f;
F = @(u) D2 + 2*diag((I*u))*I;
A = @(u) [F(u); diffrow(n+2,0,-1); diffrow(n+2,0,1)];
rhs = @(u) [N(u); 0; 0];
u = alpha*(1-x)/2 + beta*(1+x)/2;
for k = 1:8
  disp(norm(N(u))), u = u - A(u)\rhs(u); 
end
plot(chebfun(u),'.-',LW,.7,MS,16)
xlabel('$x$',FS,16,IN,LT), ylabel('$u$',FS,16,IN,LT)
set(gca,FS,11)
plottitle(9,fs,'left'), print -depsc ex09

% ex10.m - EXAMPLE 10
abbreviations
n = 80; mu = 2;
X = [0,1]; I = eye(n);
D = diffmat([n n],X,'trig');
D2 = diffmat([n n],2,X,'trig');
N = @(u,T) D2*u - T*mu*(1-u.^2).*(D*u) + T^2*u;
dNdu = @(u,T) D2 - T*mu*diag(1-u.^2)*D + 2*T*mu*diag((u).*(D*u)) + T^2*I;
dNdT = @(u,T) -mu*(1-u.^2).*(D*u) + 2*T*u;
A = @(u,T) [dNdu(u,T) dNdT(u,T); diffrow(n,0,0,X,'trig') 0];
rhs = @(u,T) [N(u,T); 0];
t = trigpts(n,X);
u = 2*cos(2*pi*(t+1/6));
T = 5;
for k = 1:9
  disp(norm(N(u,T)))
  vs = A(u,T)\rhs(u,T); 
  u = u - vs(1:n);
  T = T - vs(n+1); disp(T)
end
t = T*t; X = T*X;  % rescale to [0,T]
plot(chebfun(u,X,'trig'),'.-',LW,.7)
xlabel('$\tau$',FS,16,IN,LT), ylabel('$u$',FS,16,IN,LT)
axis([0 T -2.5 2.5]), set(gca,FS,11,'xtick',T*(0:4)/4)
plottitle(10,fs,'left'), print -depsc ex10

% ex11.m - EXAMPLE 11
abbreviations
n = 22;
Z = zeros(n,n+2); z = zeros(1,n+2);
D2 = diffmat([n n+2],2);
D = diffmat([n n+2],1);
I = diffmat([n n+2],0);
lefteval = diffrow(n+2,0,-1);
righteval = diffrow(n+2,0,1);
ep = 0.05;
A = [ Z          ep*D2+D
      ep*D2-D    Z
      lefteval   z
      righteval  z
      z          lefteval
      z          righteval ];
B = [I Z; Z I; z z; z z; z z; z z];
[UV,Lambda] = eig(A,B); d = diag(Lambda);
ii = find(d>-1e-12); 
d = d(ii); UV = UV(:,ii);
[d,ii] = sort(d);
UV = UV(:,ii); UV = UV(:,2:end); sig = d(2:end);
U = UV(1:n+2,:); V = UV(n+3:end,:);
format long
subplot(2,2,2)
for j = 1:3
  u = U(:,j); if u(end-2)<0, u = -u; V(:,j) = -V(:,j); end
  plot(chebfun(u),'.-',LW,.7,MS,7), hold on
  disp(sig(j)), ylim([-1.3 1.6])
end
xlabel('$x$',FS,12,IN,LT), ylabel('$u$',FS,12,IN,LT), set(gca,FS,8)
subplot(2,2,1)
for j = 1:3
  v = V(:,j);
  plot(chebfun(v),'.-',LW,.7,MS,7), hold on
  ylim([-1.3 1.6])
end
xlabel('$x$',FS,12,IN,LT), ylabel('$v$',FS,12,IN,LT), set(gca,FS,8)
plottitle(11,fs), print -depsc ex11

% ex12.m - EXAMPLE 12
abbreviations
n = 45; d = 10; a = 1; b = 1; alpha = 1; beta = 0;
X = [0 d];
I = diffmat([n n+2],0,X);
D2 = diffmat([n n+2],2,X);
N = @(u,w) D2*u + sin(I*u) - I*w;
M = @(u,w) D2*w + cos(I*u).*(I*w) + (a/b)*(I*u);
dNdu = @(u,w) D2 + diag(cos(I*u))*I;
dNdw = @(u,w) -I;
dMdu = @(u,w) (a/b)*I - diag(sin(I*u).*(I*w))*I;
dMdw = dNdu;
z = zeros(1,n+2);
A = @(u,w) [dNdu(u,w) dNdw(u,w)
            dMdu(u,w) dMdw(u,w)
            diffrow(n+2,0,0,X) z
            diffrow(n+2,1,0,X) z
            z diffrow(n+2,0,d,X)
            z diffrow(n+2,1,d,X) ];
rhs = @(u,w) [N(u,w); M(u,w); 0; 0; 0; 0];
t = chebpts(n+2,X); u = 1 + 0*t; w = 0*t;
for k = 1:6
  vq = A(u,w)\rhs(u,w); 
  u = u - vq(1:n+2);
  w = w - vq(n+3:end);
end
plot(chebfun(u,X),'.-',LW,.7,MS,12), hold on
plot(chebfun(w,X),'.--r',LW,.7,MS,12)
xlabel('$t$',FS,16,IN,LT), ylabel('$u,w$',FS,16,IN,LT)
set(gca,FS,11), ylim([-0.6 1.2]), grid on, hold off
plottitle(12,fs), print -depsc ex12