function p5
  a = random_pentadiag(6)
  [l,u] = crout_pentadiag(a)
  a_minus_lu = a - l*u
end
function [l, u] = crout_pentadiag(a)
  [n, dummy] = size(a);

  l = zeros(n,n);
  u = zeros(n,n);

  for j = 1:n
    for i = max(j-2,1):j-1
      u(i,j) = a(i,j);
      for k = max(1,j-2):i-1
        u(i,j) = u(i,j) - l(i,k)*u(k,j);
      end
      u(i,j) = u(i,j)/l(i,i);
    end
    u(j,j)=1;
    for i = j:min(n,j+2)
      l(i,j) = a(i,j);
      for k = max(1,i-2):j-1
        l(i,j) = l(i,j) - l(i,k)*u(k,j);
      end
    end
  end
  assert(abs(a-l*u) < 1e-15)
end
function a = random_pentadiag(n)
  a = zeros(n,n);

  for i = 1:n
    for j = max(i-2,1):min(n,i+2)
      a(i,j) = rand;
    end
  end
end
function assert(v)
  if (v == 0)
    error('Assertion failed.')
  end
end
