There is no way because there is no solution to the given system with $x_2 \neq 0$. This is because the second block of the equation system reads $Ax_2 = 0$, which has no non-trivial solutions because $A$ is SPD (Otherwise $x_2$ would be an eigenvector of $A$ corresponding to a zero eigenvalue).
Sagemath looks fast enough. Here is an example created by taking a 999x1000 random rational matrix and appending its column sums at the bottom, so that [1, 1, ... , 1, -1] is the kernel.
sage: A = random_matrix(QQ, 999, 1000)
sage: B = A.stack(sum(A))
sage: %time v = B.kernel()
CPU times: user 5.21 s, sys: 36.4 ms, total: 5.24 s
Wall time: 5.24 s
I have used the method you described for a similar problem. I found that a combination of your Newton iteration, along with Brent's method on minimizing the smallest singular value gave pretty good results. I believe I used Newton's method first for a few iterations, then switched over to Brent's. This paper of mine contains the description for the analogous ...
The problem of finding singular points of a parameter-dependent matrix is called the nonlinear eigenvalue problem in the field of numerical linear algebra. I think you could benefit from using developed methods. It is a quite active area of research. You can consider solving it with the Julia software package NEP-PACK (I am a lead developer of NEP-PACK). NEP-...
In MATLAB, null([A; B]) will find an orthogonal basis for the intersection of the null spaces of A and B.
It seems unlikely that you really want to find this basis, but it's not clear from your question what you're actually trying to do.
If your sparse matrix isn't too large and you can store it in memory, you could use a sparse (rank-revealing) QR (or LU, or SVD) factorization to determine the kernel of your matrix $A$.
Alternately, if you know something about the spectrum of your operator, you could apply a shift to it -- solve for eigenvectors associated with an eigenvalue of $\lambda$ ...