Loss function returns nan on time series dataset using tensorflow
This is because the the input values are not normalised to a standard scale. X consists of values ranging from 0 to 9999. This would make the training unstable. If the idea is to train the net to predict the sequence that you described, an easy fix is to feed it smaller values in the training data. For instance, feed it the input, 

output = split_sequence(range(10),n_steps)


Along with a batch size of say, 5. It will learn to generalize to higher numbers as well. 

With the changes made and the following for loop, the output is as expected. 

for epoch in range(n_epochs):
    for X_batch, y_batch in generate_batch(input, output, batch_size):
        # print(X_batch,y_batch)
        mse_val, _ = session.run([mse, training_op], feed_dict={
                                 X: X_batch, y: y_batch})

    if epoch % 100 == 0:
        print('epoch', epoch, 'MSE=', mse_val)
        error_list.append(mse_val)

X_test = np.array([[45, 46, 47], [57, 58, 59], [1000, 1001, 1002]])
y_test = session.run(y_predictions, feed_dict={X:X_test})

print(y_test) # Should be 48, 60, 1003

plt.plot(error_list)
plt.show()


Final y_test is 

[[  47.997036]
 [  59.996166]
 [1002.9279  ]]


If not this change, another change that you can make is to use the Adam Optimiser instead of SGD. It will make the learning more stable and it works even with bigger values.