Using LSTM's on Multivariate Input AND Multivariate Output
You can simply change the number of units at the last layer if you want to predict multiple outputs (or just said differently, one output with multiple features). For example from the link you shared, you can change the units of the Dense layer:

model = Sequential()
model.add(LSTM(50, input_shape=(train_X.shape[1], train_X.shape[2])))
model.add(Dense(2))
model.compile(loss='mae', optimizer='adam')


The above model would now predict the next step of an output with 2 "features". 
Note that you output should be of shape num_samplesx2 now.