clear;

N = 33;
h = 1 / (N - 1);
n = 0;

for r = [1/6, 1/2, 65/128]
  k = r * h ^ 2 ;
  n = n + 1;
  M = round(300 / r);
  u = zeros(N,M);
  u(N,1) = 1;
  P = toeplitz([1 - 2 * r, r, zeros(1, N - 2)]);
  for m = 1:(M - 1)
    u(:,m + 1) = P * u (:,m);
    u(1,m + 1) = 0;
    u(N,m + 1) = 1;
  end
  dk = round(20 / r);
  uu = u(:, 1:dk:M);
  [X Y] = meshgrid([0:dk * k:(M - 1) * k], [0:h:1]);
  subplot(2,2,n);
  plot(Y(:,1:2:(M / dk)), uu(:,1:2:(M / dk)), '-b');
end
subplot(2,2,4);
mesh(X,Y, uu);