How to Predict the future values of time horizon with Keras?
This answer goes a little bit in a different direction, but I hope it still answers your question. It uses the idea of a rolling forecast/prediction.

Because you use the word horizon, I will assume you mean that you would like to predict 10 days into the future at a given time step. There are a few ways of doing this. With this kind of time-series problem, it is common to make the assumption that only a certain history will influence the next few time steps (neglecting seasonal effects). 

Example in words:

So in your case, you might use e.g. the previous 60 days, and predict the next 10.
Taking your 100 rows of data as an example, this means you can actually make (100 - 60 - 9) = 31 predictions, each prediction of 10 time steps ahead (we will need these 31 predictive_blocks later). From 100 rows we lose the first 60 to fit the first model. Of the remaining 40 rows of data, we can predict 10 steps ahead (rows 61-70), then we shift the whole thing one row further and repeat. The last prediction of 10 future points would be for rows 91-100. After that we cannot predict 10 steps anymore, so we stop - and this is why we have to subtract that extra 9. [There are of course ways to continue making prediction, as to use all the data]

Example with a thousand words:

Let me paint the picture; to help explain the idea of a shifting window prediction.

For each train set (e.g. from t=0 to t=5 in red - train set 1), you want to predict the following H time steps (corresponding to t=6 in orange - test set 1). In this, your horizon is simply one i.e. H=1.



From what I understand, you would like to predict the next 10 days, meaning you need H=10.

In order to try this with your example, I think you will need to make two changes.

Change #1

The shape of your train and test sets will need to match the new horizon. Each sample of your model input (the x_train and x_test can stay the same as before. However, each sample in your test set will have to contain the next H=10 values of the label, not just a single value.

Here is a rough example of how you might do this:

# Define our horizon
H = 10

# Create data split, using values from my example above
window_size = 60
num_pred_blocks = 31    # as computed above

# Loop over the train and test samples to create the sliding window sets
x_train = []
y_train = []
for i in range(num_pred_blocks):
    x_train_block = x_train[i:(i + window_size)]    # 31 blocks of 60 * num-columns
    x_train.append(x_train_block)
    y_train_block = y_train[(i + window_size):(i + window_size + H)]    # 31 blocks of 10 * 1
    y_train.append(y_train_block)


Because you are doing out-of-sample testing, your predictions are already interesting to look analyse. Once this runs, you can then create the equivalent test datasets with the new data you mentioned.

Without knowing your data too well, I don't know if your should be predicting the y-values of the same row as the input, or of the following row. Additionally, depending on your data, you could be including the past values of y in each of the x_train blocks. In this case you'd simply swap x for the whole table i.e. data[cols], where new_cols = ['Demand'] + cols.

Change #2

You will need to make the model reflect this horizon, by forcing it to output H values.

Here is an example of how to specify the model:

# Define our horizon
H = 10

# Create the model using the parameterised horizon
fit1 = Sequential ()
fit1.add(LSTM(output_dim = 4, activation='tanh', input_shape =(4, 1)))
fit1.add(Dense(output_dim=30, activation='sigmoid')
fit1.add(Dense(output_dim=H))    # our horizon is produced!


Note:
In your model specification, you don't need to add the final linear Activation, as the preceding Dense layer by default includes a linear activation. See the excellent documentation here.

This is a big topic and there are many things that you could try out. I agree with the comments on your question, that you will need a lot more data to allow an RNN to make a meaning representation of the model.

If you are not just doing this to learn about LSTMs etc., another practical approach might be to look into simpler time-series models such as an ARIMA model (do not be intimidated by the complicated name - it is much simpler than an LSTM). Such models can be constructed quite easily with Python, using the statsmodels package, which has a nice implementation.