  function [x0,X] = MarkovSol(f,g,u,v,b)
  % f, g, and b are column n-vectors
  % u and v are column q-vectors
  % x0 = A(f,g)\b
  % X is an n-by-q matrix with the property that if e1 is the first column of eye(n,n)
  % then X(:,j) = A(f+u(j)*e1,g+v(j)*e1)\b for j=1:q
  
  n = length(f); 
  % Set up C = C(f,g) and d = d(f,g)...
  C = zeros(n,n);
  for i=1:n
      for j=i+1:n
          C(i,j) = 1/ ((f(i)-f(j))^2 + (g(i)-g(j))^2);
      end
  end
  C = C + C';
  d = C*ones(n,1);
  
  % We will "live" off of this...
  [L,U,P] = lu(C);
  
  % Fact 1. Solving systens with C and D = diag(d)
  %
  % To solve b = Ax = C*inv(D)x we first solve C*xtilde = b then set
  %      x = D*xtilde
  
  % Fact 2. Connection between C(f,g) and Ctilde = C(f+alfa*e1,g+beta*e1)
  %
  % Let s be the first column of C(f,g), i.e.,
  %      s = [0 ; 1./(f(1)-f(2:n))^2 + (g(1) - g(2:n)).^2]
  % Let t be the first column of C(f+alfa*e1,g+beta*e1), i.e.,
  %      t = [0 ; 1./(f(1)+alfa-f(2:n))^2 + (g(1)+beta - g(2:n)).^2]     
  % Note that
  %      Ctilde =  C(f+alfa*e1,g+beta*e1) 
  %             = C(f,g) - s*e1'- e1*s' + t*e1' + e1*t'
  %             = C(f,g) +(t-s)e1' + e1*(t-s)'
  %             = C + W*Z'
  % where W = [t-s e1]  and Z = [e1 t-s]. 
  
  % Fact 3. Connection between d(f,g) and dtilde = d(f+alfa*e1,g+beta*e1)
  %
  % How do the column sums change when we change f(1) and g(1)?
  % The whole first column of C changes but only the top entry changes
  % in columns 2 through n. Thus,
  %      dtilde = d(f+alfa*e1,g+beta*e1) = [sum(t); d(2:n)-s(2:n)+t(2:n)]
  
  % Fact 4. exploiting the connection between inv(C) and inv(C+W*V')
  %
  % The Sherman-Morrison-Woodbury formula says
  %       inv(C+W*Z') = inv(C) - inv(C)*W*inv(M)*Z'*inv(C)   M = I + Z'*inv(C)*W
  % so if Cr = b then the solution to (C+W*Z')q = b is given by
  %
  %       q = inv(C+WZ')b = r - inv(C)*W*inv(M)*Z'*r
  %                       = r - U\(L\(P*(W*(M\(Z'*r)))))
  % where M = I + Z'*(U\(L\(P*W)))
  
  % Solve A(f,g)*x0 = b;
  r = U\(L\(P*b));
  x0 = d.*r;
  
  m = length(u);
  X = zeros(n,m);
  s = C(:,1);
  e1 = [1;zeros(n-1,1)];
  for k=1:m
      alfa = u(k);
      beta = v(k);
      % Solve A(f+alfa*e1,beta*e2)*X(:,k) = b
      t = [0 ; 1./(   (f(1)+alfa-f(2:n)).^2   + (g(1)+beta - g(2:n)).^2  )] ;
      W = [t-s e1];
      Z = [e1 t-s];
      M = eye(2,2) + Z'*(U\(L\(P*W)));
      q = r - U\(L\(P*(W*(M\(Z'*r)))));
      dtilde = [sum(t); d(2:n)-s(2:n)+t(2:n)];
      X(:,k) = dtilde.*q;
  end
      
      