Time-series multi-step generalization from single step model
Welcome to the community Fra, below you can find a worked out example implementing a multivariate several input features (as I think is your case) time series forecasting, predicting multiple future steps (multi-step forecast), applying bayesian hyperparametrization.
It is based on a from simpler to more complex approach, so you can see there are few layers in the architecture to begin with, and what is being hyperparametrized is the units in such layers. This is a whole worked out example, in case you can almost directly apply it and try out modelling first with a single output node for a single-step prediction, and then add more future step predictions to check if you still run on the same issue as you described (beggining with n_output = 1 and then incrementing this, e.g. in my case it is 72 future values to predict).
For a dataset of this type:

and trying to predict several future time steps (the deep blue dots) for the temperature values:

we have some code like:

train-validation dataset splitter functions:
  def multivariate_data(dataset, target, start_index, end_index, history_size, target_size, step, single_step=False):
    data = []
    labels = []

    start_index = start_index + history_size
    if end_index is None:
      end_index = len(dataset) - target_size

    for i in range(start_index, end_index):
      indices = range(i-history_size, i, step)
      data.append(dataset[indices])

      if single_step:
        labels.append(target[i+target_size])
      else:
        labels.append(target[i:i+target_size])

    return np.array(data), np.array(labels)



This is what we get by applying that function (in your case, you do not need to apply the sampling interval of 6 timesteps, you can set it to 1; this is in case you want to downsample a maybe too long history):


function implementation to build the train-eval sets:
  past_history = 720 # historic values to consider for each prediction
  future_target = 72 # number of future values to predict 
  STEP = 6 # sampling frequency (it makes taking 720/6 historic points for each sample)

  x_train_multi, y_train_multi = multivariate_data(dataset_standardized, dataset_standardized[:, 1], 0,TRAIN_SPLIT, past_history, future_target, STEP) 
  x_val_multi, y_val_multi = multivariate_data(dataset_standardized, dataset_standardized[:, 1], TRAIN_SPLIT, None, past_history, future_target, STEP)



This way you have:

Single window of past history : (120, 3)
Target temperature to predict : (72,)

For the training phase, we define the batch size and other values and build our training and validation sets:
BATCH_SIZE = 256
BUFFER_SIZE = 10000

train_data_multi = tf.data.Dataset.from_tensor_slices((x_train_multi, y_train_multi))
train_data_multi = train_data_multi.cache().shuffle(BUFFER_SIZE).batch(BATCH_SIZE).repeat()

val_data_multi = tf.data.Dataset.from_tensor_slices((x_val_multi, y_val_multi))
val_data_multi = val_data_multi.batch(BATCH_SIZE).repeat()

and defining a simple architecture to begin with, we have:
multi_step_model = tf.keras.models.Sequential()
multi_step_model.add(tf.keras.layers.LSTM(32,
                                          return_sequences=True,  # &quot;Boolean. Whether to return the last output. in the output sequence, or the full sequence&quot;
                                          input_shape=x_train_multi.shape[-2:]))
multi_step_model.add(tf.keras.layers.LSTM(16, activation='relu'))
multi_step_model.add(tf.keras.layers.Dense(72))

multi_step_model.compile(optimizer=tf.keras.optimizers.RMSprop(clipvalue=1.0), loss='mae')

for which we can check the performance on the validation set as:
EVALUATION_INTERVAL = 200
EPOCHS = 10

multi_step_history = multi_step_model.fit(train_data_multi, epochs=EPOCHS,
                                          steps_per_epoch=EVALUATION_INTERVAL,
                                          validation_data=val_data_multi,
                                          validation_steps=50) 

having the predictions in red below, after de-scaling the values to its original scale:
for x, y in val_data_multi.take(3):
    x_denormalized = (x[0])*temps_std+temps_mean
    y_denormalized = (y[0])*temps_std+temps_mean
    
    multi_step_plot(x_denormalized, y_denormalized, (multi_step_model.predict(x)[0])*temps_std+temps_mean)


and maybe the most interesting part for your use case, you can try implementing hyperparametrization of some params as hidden layers units as follows with the keras bayesian tuner:
from kerastuner import BayesianOptimization

def build_model_2(hp):
    model = keras.Sequential()
    model.add(keras.layers.LSTM(units=hp.Int('units',min_value=16,
                                        max_value=32,
                                        step=16), 
                                        #activation='relu', 
                                        return_sequences=True, 
                                        input_shape=x_train_multi.shape[-2:]))
    model.add(keras.layers.LSTM(units=hp.Int('units',min_value=16,
                                        max_value=32,
                                        step=16),  
                                        activation='relu'))
    model.add(keras.layers.Dense(72))

    model.compile(loss='mae', optimizer=tf.keras.optimizers.RMSprop(clipvalue=1.0),
            
                   metrics=['mae'])
    return model

# define model
bayesian_opt_tuner_2 = BayesianOptimization(
    build_model_2,
    objective='mae',
    max_trials=3,
    executions_per_trial=1,
    directory=os.path.normpath('C:/keras_tuning'),
    project_name='timeseries_temp_ts_test_from_TF_ex_multivar_multistep_2',
    overwrite=True)

EVALUATION_INTERVAL = 200
EPOCHS = 10
bayesian_opt_tuner_2.search(train_data_multi, 
             epochs=EPOCHS,
             validation_data=val_data_multi,
             validation_steps=50,  
             steps_per_epoch=EVALUATION_INTERVAL)

# and selecting the best model: 

    best_MULTIVAR_MULTISTEP_LSTM_model_2 = bayesian_opt_tuner_2.get_best_models(num_models=1)[0]

we have our new model predictions:

which gave an improved MAE on the validation set.