The first thing you need to remember is that when you ask to solve() a vector of equations, then solve() has to find a solution that satisfies all of the entries simultaneously. Because your x is length 1000 and there are two entries for each x, you are asking to solve 2000 simultaneously equations.
Secondly, if you use a symbolic x, and go through doing stepwise elimination of the variables, you will find that indeed 0 are the only values that will reliably make that pair of equations 0 on both sides. That does not mean there are no other zeros, but the other zeros are in x, not in b or tau_new . The additional zero is at x = exp(-LambertW(20000*sqrt(2)*Pi*exp(1/3)*(1/5421))-2/3) which is about 0.05157587627 which is before the start of your x search range.