Allowing for second feature in BiLSTM Model

I am working on a python code for stock market prediction. I am hoping to update the BiLSTM model to allow for 2 or more features/columns automaically in the inputs to improve model performance. (Eg: Prophet model yhat results.)

I would appreciate any help in trying to figure this out as I am unsure on where i should update to allow for this. Would any additional columns also be scaled for the model to work most effieciently?

I have tried to update the below but am getting an error regarding the operand combination

Data:

Decomp'd BiLSTM

def de_BiLSTM(train_dataset, test_dataset, n_input = 3, n_features = 2, u = 32, epochs = 100, batch = 1, verbose = 0):
    try:
        if (verbose == 1):
            print("starting BiLSTM")
        train=pd.DataFrame(data=train_dataset["y"],columns=["y"])
        if (verbose == 1):
            print("starting BiLSTM2")
        start = train_dataset.ds[0]
        if (verbose == 1):
            print("start value:", start)
        dti = pd.date_range(start, periods=len(train_dataset), freq='D')
        train["time"]=dti
        train_days=pd.DataFrame()
        train_days["time"]=dti
        train.set_index('time', inplace=True)
        if (verbose == 1):
            print("starting testdataset")
        test=pd.DataFrame(data=test_dataset["y"],columns=["y"])
        if (verbose == 1):
            print("starting testdataset2", test_dataset.head(1))
        #start2 = test_dataset.ds[0]
        start2 = test_dataset['ds'].iloc[0]
        if (verbose == 1):
            print("start value:", start2)
        dti2 = pd.date_range(start2, periods=len(test_dataset), freq='D')
        test["time"]=dti2
        if (verbose == 1):
            print("test_days")
        test_days=pd.DataFrame()
        test_days["time"]=dti2
        test.set_index('time', inplace=True)

        if (verbose == 1):
            print("starting decomposition")
        res = sm.tsa.seasonal_decompose(x=train["y"], extrapolate_trend='freq')
        lstm_train = pd.DataFrame()
        lstm_train["y"]=res.trend+res.resid

        res_test = sm.tsa.seasonal_decompose(x=test["y"], extrapolate_trend='freq')
        lstm_test = pd.DataFrame()
        lstm_test["y"]=res_test.trend+res_test.resid
        
        if (verbose == 1):
            print("starting scaler")
        scaler = MinMaxScaler()
        scaler.fit(lstm_train)
        scaled_train_data = scaler.transform(lstm_train)
        scaled_test_data = scaler.transform(lstm_test)
        
        if (verbose == 1):
            print("starting generator")
        generator = TimeseriesGenerator(scaled_train_data, scaled_train_data, length=n_input, batch_size=batch)
        lstm_model = models.Sequential()
        lstm_model.add(layers.Bidirectional(layers.LSTM(units=u, return_sequences = True, activation='relu', input_shape=(n_input, n_features))))
        lstm_model.add(layers.Dropout(0.005))
        lstm_model.add(layers.Bidirectional(layers.LSTM(units=u, return_sequences = True)))
        lstm_model.add(layers.Dropout(0.005))
        lstm_model.add(layers.Bidirectional(layers.LSTM(units=u, return_sequences = True)))
        lstm_model.add(layers.Dropout(0.005))
        lstm_model.add(layers.Bidirectional(layers.LSTM(units=u)))
        lstm_model.add(layers.Dropout(0.005))
        lstm_model.add(layers.Dense(1))
        lstm_model.compile(optimizer='adam', loss='mse')
        if (verbose == 1):
            print("starting fitting")
        lstm_model.fit_generator(generator,epochs=epochs, verbose=verbose)
        batch = scaled_train_data[-n_input:]
        current_batch = batch.reshape((1, n_input, n_features))
        lstm_predictions_scaled = list()

        for i in range(len(lstm_test)):   
            lstm_pred = lstm_model.predict(current_batch)[0]
            lstm_predictions_scaled.append(lstm_pred) 
            current_batch = np.append(current_batch[:,1:,:],[[lstm_pred]],axis=1)

        lstm_predictions = scaler.inverse_transform(lstm_predictions_scaled)
        lstm_result = pd.DataFrame()
        lstm_result['y'] = test_dataset['y']
        lstm_result['LSTM_Predictions'] = lstm_predictions
        
        
        l_r2 = r2_score(lstm_result.y, lstm_result.LSTM_Predictions)
        l_mse = mean_squared_error(lstm_result.y, lstm_result.LSTM_Predictions)
        l_mae = mean_absolute_error(lstm_result.y, lstm_result.LSTM_Predictions)
        l_rmse = math.sqrt(l_mse)
        if (verbose==1):
            print("r2:", l_r2, "mse:", l_mse, "rmse:", l_rmse, "mae:", l_mae)
        
        
        seasonal_lstm = pd.DataFrame()
        seasonal_lstm["y"] = lstm_train["y"] + res.seasonal

        seasonal_lstm_test = pd.DataFrame()
        seasonal_lstm_test["y"] = lstm_test["y"] + res_test.seasonal

        seasonal_lstm_pred = pd.DataFrame()
        #seasonal_lstm_pred["y"] = lstm_result['LSTM_Predictions'] + res_test.seasonal
        seasonal_lstm_pred["s"] = res_test.seasonal  
        seasonal_lstm_pred["p"] = [x for x in lstm_result['LSTM_Predictions']]
        seasonal_lstm_pred["y"] = seasonal_lstm_pred[['s', 'p']].sum(axis=1)
        seasonal_lstm_pred.drop('s', inplace=True, axis=1)
        seasonal_lstm_pred.drop('p', inplace=True, axis=1)

        #seasonal_lstm_pred["y"] = seasonal_lstm_pred["y"] + lstm_result['LSTM_Predictions']
        return lstm_result,seasonal_lstm, seasonal_lstm_test, seasonal_lstm_pred, l_r2, l_mae, l_mse, l_rmse
        #return lstm_model, lstm_train, lstm_test, lstm_result, train_days, test_days, current_batch, scaler
    except Exception as e:
        print("exception:", e)
        
   

'''