%function PlotSolutionSpace(A,b)

	close all
	clear
	figure; hold on; 

A = ones(3);
A = rand(3); A(:,3) = A(:,1)+A(:,2)

% sample point on span of col(A)
u = A*[1; 1; 1];

% sample point on span of col(A)
v = A*[1; 2; 3];

% special point on span of col(A), which
% we'd like to know the linear combination
% of col vectors of A necessary to find this vectorA


c = 2*(-1+2*rand(3,1)); 
b = A*c;
%b = 1/norm(b) * b;
b = [1 1 0]'
b0 = b


M = [A b];
R = rref(M);

	% The origin
	plot3(0,0,0,'.k','MarkerSize',20);
	h0 = plot3(0,0,0,'ok','MarkerSize',20);
	
	plot3(b(1),b(2),b(3),'.k','MarkerSize',20);
	hb = plot3(b(1),b(2),b(3),'sk','MarkerSize',20);

	pause
if rank(A) == 1

	% A is rank 1, col(A) = 1-D
	a = sum(A,2);
	imageA = a*[-10:.01:10];
	hs = plot3(imageA(1,:),imageA(2,:),imageA(3,:),'.k');

	disp('The column span of A')
	pause

	% A is Rank 1, null space has dim 2
	syms y z
	x0 = -R(1,2)*y - R(1,3)*z
	ezmesh(x0,[-10,10,-10,10]);

	for i = 1:100;
			X0(:,i) = subs([x0;y;z],{y,z},{-10+20*rand(1),-10+20*rand(1)});
	end
	
	disp('The solution to Ax=0')
	pause
	% Soln space to Ax = b

	xb = -R(1,2)*y - R(1,3)*z + R(1,4)
	ezmesh(xb,[-10,10,-10,10]);
	disp('The solution to Ax=b')


	% Check accuracy
	for i = 1:100;
			Xb(:,i) = subs([xb;y;z],{y,z},{-10+20*rand(1),-10+20*rand(1)});
			N(i) = norm(A*[subs([xb;y;z],{y,z},{-10+20*rand(1),-10+20*rand(1)})] - b);
			if N>0.1; disp('Error in xb'); return; end;
	end

	

	%legend([h0,hb,hs,hX],'The origin','The vector b','Col ( A )', 'Soln to Ax=b','Soln to Ax=0','Location','BestOutside')

elseif rank(A) == 2


	% A is Rank 1, soln to Ax = b has dim 2 
	
	
	% A is Rank 2, soln to Ax = b has dim 1
	x_soln = [-R(1,3); -R(2,3); 1];
	
	X = x_soln*[0:.01:10];
	
	% the solution space to Ax = b is the line
	% through x_b where Ax_b = b, parallel to X below
	% Xb = b + X

	syms z
	y = -R(2,3)*z + R(2,4);
	x = -R(1,2)*y - R(1,3)*z + R(1,4);

	xb = subs([x;y;z]',{z},1)';
		
	Xb = repmat(xb,1,size(X,2)) + X;
	
	
	
	%PlotPlane(b,u,v);
	
	PlotPlane;
	
	pause
	
	hXb = plot3(Xb(1,:),Xb(2,:),Xb(3,:),'.m','MarkerSize',20);
	
	
	hX = plot3(X(1,:),X(2,:),X(3,:),'.k','MarkerSize',20);
	legend([h0,hb,hXb,hX],'The origin','The vector b','Soln to Ax=b','Soln to Ax=0')
end

b = Xb(:,end)
b0 
plot3(b(1),b(2),b(3),'db','MarkerSize',20);
for i = 1:500;

	plot3(b(1),b(2),b(3),'db','MarkerSize',20);
	hb = plot3(b(1),b(2),b(3),'sk','MarkerSize',20);
	drawnow;
	pause(.1)
	b = A*b;
end
