clear;

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

for r = [1 2 4]
  k = r * h ^ 2;
  n = n + 1;
  M = round(4800 / r);
  
  u = zeros(N,M);
  u(N,1) = 1;
  
  C = toeplitz([-2, 1, zeros(1,N - 2)]);
  A = sparse(2 * eye(N,N) - r * C);
  B = sparse(2 * eye(N,N) + r * C);
  b = zeros(N,1);
  b($) = 2 * r;
  
  for m = 1:(M - 1)
    u(:, m + 1) = A \ (B * u(:,m) + b);
    u(1, m + 1) = 0;
    u(N, m + 1) = 1;
  end
  
  dk = 80;
  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);