function x_n = qsd(A,b)
	Asz = size(A);
	bsz = size(b);
	assert(Asz(1) == Asz(2) && Asz(1) == bsz(1));
	% A should also be symmetric
	x = 0;
	while true
		a = (b'*b)/(b'*A*b);
		x_n = x + a*b
        if abs(x_n - x ) < .0001
			break;
		else
			x = x_n;
		end
		b = b - a*(A*b);
	end
end
	
