Same probability for all classes
Assuming that you normalized the pixel-values as it is done in the tutorial, your inputs were vectors of numbers betweeen 0 and 1. Now, if your weight matrices are also randomly drawn numbers between 0 and 1, the input to the hidden layer will be sums of 783 numbers between 0 and 1, i.e. probably something > 100. Now, check out the sigmoid function and its derivative



As you see, it saturates quite quickly for input values > 5. If you choose the initialization as you did, all hidden neurons should be very close to one, while at the same time the derivative should be very close to zero. This would explain how all outputs of the softmax were equally close to 0.1 and since the gradient was close to zero, the network probably didn't learn anything.

Once you changed the weight initialization to something between -1 and 1, the inputs to the hidden layer should be sums of numbers fairly evenly distributed around 0, thus the sigmoid output was about 0.5 and, most importantly, the gradient was non-zero so your network actually got trained.

As you already noticed, choosing the initialization wisely is crucial for getting decent results. Initializations should also take the number of input neurons into account, otherwise you initialize the network in such a way, that the gradients will be close to zero.

A similar problem can occur if you do not normalize the input data properly, e.g. if you feed the pixel-values between 0 and 255 directly to the model.

I'm not sure about papers on this, but maybe start with the original work on Glorot initialization or checkout the initializers Tensorflow has to offer. On that site under "Functions" they list common initializers and also link to the respective papers.

I hope this was helpful.