  N = 40;
  mu = 1.0e-3;
  h  = 1/N;
  N1    = N+1;
  N1N1  = N1*N1;
  N4    = round(N/4);
  ustar = [zeros(N4,1);ones(N1-2*N4,1);zeros(N4,1)];
  ustar = kron(ustar,ustar');
  v     = ustar + 0.1*randn(N1,N1);

  v = reshape(v,N1N1,1);

  Dxdiag = -ones(N1,N1); Dxdiag(N1,:) = 0;
  Dxsup  = ones(N1,N1);  Dxsup(1,:)   = 0;
  Dx = spdiags(Dxdiag(:),0,N1N1,N1N1) ...
     + spdiags(Dxsup(:),1,N1N1,N1N1);
  Dydiag = -ones(N1,N1); Dydiag(:,N1) = 0;
  Dysup  = ones(N1,N1);  Dysup(:,1)   = 0;
  Dy = spdiags(Dydiag(:),0,N1N1,N1N1) ...
     + spdiags(Dysup(:),N1,N1N1,N1N1);
  D = Dx'*Dx + Dy'*Dy;
  A = mu*D/h^2 + speye(N1N1);

tic
  u = A \ v;
toc

  u = reshape(u,N1,N1);
  v = reshape(v,N1,N1);

u1 = u;

  h1 = figure(1); clf; set(h1,'Position',[10 10 900 300]);
  subplot(1,3,1); 
  imagesc(ustar); colormap('gray'); axis image; axis off;
  title('exakt')
  subplot(1,3,2); 
  imagesc(v); colormap('gray'); axis image; axis off;
  title('verrauscht')
  subplot(1,3,3); 
  imagesc(u); colormap('gray'); axis image; axis off;
  title('entrauscht')

tic
  v = reshape(v,N1N1,1);
  Av = full([A,v]);
  n  = N1N1;
  for k=1:(n-1)
    mm = min(n,k+N1);
    for i=(k+1):mm
      mp = Av(i,k)/Av(k,k);
      for j=(k+1):mm
        Av(i,j) = Av(i,j) - mp*Av(k,j);
      end
      Av(i,n+1) = Av(i,n+1) - mp*Av(k,n+1);
    end
  end

  u = zeros(n,1);
  u(n) = Av(n,n+1)/Av(n,n);
  for i=(n-1):(-1):1
    mm = min(n,i+N1);
    s = sum(Av(i,(i+1):mm)*u((i+1):mm));
    u(i) = (Av(i,n+1)-s)/Av(i,i);
  end
  u = reshape(u,N1,N1);
  v = reshape(v,N1,N1);
  disp(strcat('diff in solutions=',num2str(norm(u(:)-u1(:)))))
toc

tic
  v = reshape(v,N1N1,1);
  Av = full([A,v]);
  n  = N1N1;
  for k=1:(n-1)
    mm=min(n,k+N1);
    for i=(k+1):mm
      Av(i,k) = Av(i,k)/Av(k,k);
      Av(i,(k+1):mm) = Av(i,(k+1):mm) - Av(i,k)*Av(k,(k+1):mm);
      Av(i,n+1) = Av(i,n+1) - Av(i,k)*Av(k,n+1);
    end
  end

  u = zeros(n,1);
  u(n) = Av(n,n+1)/Av(n,n);
  for i=(n-1):(-1):1
    mm = min(n,i+N1);
    s = sum(Av(i,(i+1):mm)*u((i+1):mm));
    u(i) = (Av(i,n+1)-s)/Av(i,i);
  end
  u = reshape(u,N1,N1);
  v = reshape(v,N1,N1);
  disp(strcat('diff in solutions=',num2str(norm(u(:)-u1(:)))))
toc
