I hope I got all the numbers into Mathematica correctly. I calculated a least squares solution which is shown in the first image. You might check it against what Matlab gives. Also give a quick look over at the data in my image.
To solve this system with constraints, I used a command "NMinimize" in Mathematica. It can minimize nonlinear functions.
In the second image you can see in the NMinimize command that I have an expression for the norm of the residual. I calculated it as the matrix product of A times a vector of Xi minus the vector b, then took the Frobenius norm. There's a variable named Bound and you can see the 10 constraints of the Xi to be between -Bound and +Bound. The result of this command is a list in curly braces {}. The first number is the norm of the residual, followed by a list of the Xi values.
Looking at the least squares solution from the first image, we see that the largest Xi, namely X4, is 19007.1. One would think that if we allow NMinimize to search for a minimum with Bound set to slightly larger than X4, it should find a solution (essentially) equal to the least squares solution. This is what it did.
In the third image, I had NMinimize find solutions for several values of Bound. I also used the plot function to plot the value of the residual norm for a range of values of Bound from 300 to 20000. It gets a little noisy as Bound approaches 20000 because I limited the number of iterations in the search so I could get some results today instead of waiting till next week. There's a lot of number crunching involved in that plot!
As you can see, the attainable residual norm when you constrain Bound to be small is substantially larger than the minimum of the least squares solution.
You posted this in the homework help forum; is this actually a homework problem? It looks like a curve fitting problem to me. I calculated the condition number of your A matrix. It's shown in the second image as the ratio of the largest to the smallest singular value, a fairly large condition number of about 36000, even though the matrix itself has a fairly small norm.
If this is a curve fitting problem, or perhaps a multiple regression problem, you may get better results with a reduced model.