  function [x1,x2] = Under2(A1,b1,A2,b2)
 % A1 is m1-by-n with rank(A1) = m1 < n, b1 is m1-by-1
 % A2 is m2-by-n with rank(A2) = m2 < n, b2 is m2-by-1
 % A1*x1 = b1 and A2*x2 = b2 and no other pair of solutions to these
 % two underdetermined systems are closer to each other in 2-norm. That
 % is, if A1*w1 = b1 and A2*w2 = b2, then norm(x1-x2,2) <= norm(w1-w2,2).
 
 % Obtain particular solution and
 % a basis for the matrix nullspace
 
 [m1,n] = size(A1);
 [Q1,R1] = qr(A1'); 
 % b = A1*x1 = A1*Q1*Q1'*x = [R1' 0][y1;z1] = R1'*y1
 y1 = R1(1:m1,1:m1)'\b1;
 x1p = Q1(:,1:m1)*y1; 
 N1 = Q1(:,m1+1:n);
 % All solutions have the form x1p + N1*z1 where z1 is (n1-m1)-by-1
 
 % Likewise...
 [m2,n] = size(A2);
 [Q2,R2] = qr(A2'); 
 y2 = R2(1:m2,1:m2)'\b2; 
 x2p = Q2(:,1:m2)*y2; 
 N2 = Q2(:,m2+1:n);
 
 
 % Find z1 and z2 to minimize 
 %                              norm( (x1p+N1*z1)-(x2p+N2*z2), 2)
 %                           =  norm( [N1 -N2]*[z1;z2] - (x2p-x1p) )
 z = [N1 -N2]\(x2p-x1p);
 z1 = z(1:n-m1);
 z2 = z(n-m1+1:2*n-m1-m2);
 x1 = x1p + N1*z1;
 x2 = x2p + N2*z2;
 
 
 
 