problem in multi step dense time series prediction using LSTM network

Hi everyone

For some time now, I have been trying to predict time series with multi features for input using the LSTM neural network. (multi input - one output)

The neural network is well trained and the error of training and test data is low. (About 1.5%) But when I want to predict the test data one by one, the error is very high and the end result is very wrong.

I put my code at the bottom. Please guide me in this regard.

clc
clear

%% Load Data
load('C:\Program Files\Polyspace\R2019b\bin\ontario_demand');
data=ontario_demand;


%% Create Train, Validation and test data
nData = numel(data);
TrainPercent = 0.995;
TestPercent = 0.005;

nTrain = round(TrainPercent*nData);
nTest = nData - (nTrain);

TrainData = data(1:nTrain,:);
TestData = data((nTrain + 1):end,:);


%% Standardize Data
mu = mean(TrainData);
sigma = std(TrainData);

StandardTrainData = (TrainData-mu)/sigma;
StandardTestData = (TestData-mu)/sigma;


%% Create X & Y for Train and Test Data
Lags = 1:5;
[XTrain, YTrain] = iLag(StandardTrainData, Lags);
[XTest, YTest] = iLag(StandardTestData, Lags);


%% Convert Matrix to Cell
XcTrain = iMakeCell(XTrain);
YcTrain = iMakeCell(YTrain);

XcTest = iMakeCell(XTest);
YcTest = iMakeCell(YTest);


%% Define LSTM Network Architecture
numFeatures = size(XTrain,1);
numResponses = size(YTrain,1);
numHiddenUnits = 100;

layers = [ ...
    sequenceInputLayer(numFeatures)
    lstmLayer(numHiddenUnits)
    fullyConnectedLayer(numResponses)
    regressionLayer];

options = trainingOptions('adam', ...
    'MaxEpochs',100, ...
    'GradientThreshold',1, ...
    'InitialLearnRate',0.05, ...
    'LearnRateSchedule','piecewise', ...
    'LearnRateDropPeriod',50, ...
    'LearnRateDropFactor',0.2, ...
    'Verbose',1, ...
    'VerboseFrequency',10, ...
    'Plots','training-progress');


%% Train LSTM Network
net = trainNetwork(XcTrain,YcTrain,layers,options);


%% Calculate Training set error indicators
[net,YPredTrain] = predictAndUpdateState(net,XcTrain);

YPredTrain = cell2mat(YPredTrain);
YPredTrain = YPredTrain*sigma + mu;
YPredTrain = YPredTrain';

YTrain2 = YTrain*sigma + mu;

RMSE_Train = sqrt(mean((YPredTrain-YTrain2).^2));
MAE_Train = mean(abs(YPredTrain-YTrain2));
MAPE_Train = mean(abs((YPredTrain-YTrain2)./YTrain2))*100;

disp('Training set error indicators :')
disp(['RMSE = ', num2str(RMSE_Train)])
disp(['MAE = ', num2str(MAE_Train)])
disp(['MAPE = ', num2str(MAPE_Train)])
disp('===================================')


%% Calculate Test set error indicators
net = resetState(net);
[net,YPredTest] = predictAndUpdateState(net,XcTest);
YPredTest = cell2mat(YPredTest);
YPredTest = sigma*YPredTest + mu;
YPredTest = YPredTest';

YTest2 = YTest*sigma + mu;

RMSE_Test = sqrt(mean((YPredTest-YTest2).^2));
MAE_Test = mean(abs(YPredTest-YTest2));
MAPE_Test = mean(abs((YPredTest-YTest2)./YTest2))*100;

disp('Test set error indicators :')
disp(['RMSE = ', num2str(RMSE_Test)])
disp(['MAE = ', num2str(MAE_Test)])
disp(['MAPE = ', num2str(MAPE_Test)])
disp('===================================')


%% Predict Test Data
net = resetState(net);

[net,YPred(1)] = predictAndUpdateState(net,XTest(:,1));
[net,YPred(2)] = predictAndUpdateState(net,[YPred(1);XTest(2:end,2)]);
[net,YPred(3)] = predictAndUpdateState(net,[YPred(2);YPred(1);XTest(3:end,3)]);
[net,YPred(4)] = predictAndUpdateState(net,[YPred(3);YPred(2);YPred(1);XTest(4:end,4)]);
[net,YPred(5)] = predictAndUpdateState(net,[YPred(4);YPred(3);YPred(2);YPred(1);XTest(end,5)]);

for i=6:size(XTest,2)
    [net,YPred(i)] = predictAndUpdateState(net,[YPred(i-1);YPred(i-2);YPred(i-3);YPred(i-4);YPred(i-5)]);
end

YPred = YPred*sigma + mu;

RMSE_Pred = sqrt(mean((YPred-YTest2).^2));
MAE_Pred = mean(abs(YPred-YTest2));
MAPE_Pred = mean(abs((YPred-YTest2)./YTest2))*100;

disp('Test set error indicators :')
disp(['RMSE = ', num2str(RMSE_Pred)])
disp(['MAE = ', num2str(MAE_Pred)])
disp(['MAPE = ', num2str(MAPE_Pred)])
disp('===================================')


%% Plot the results
figure()
plot(YTest2,'Color',[0.4 0.4 0.4])
hold on
plot(YPredTest,'Color',[0.6350 0.0780 0.1840],'LineStyle','--','LineWidth',1.5)
legend('YTest','YOutput','Location','northwest')
title(['MAE = ',num2str(MAE_Test),' , RMSE = ',num2str(RMSE_Test),' , MAPE = ',num2str(MAPE_Test),' %'])

figure()
plot(YTest2,'Color',[0.4 0.4 0.4])
hold on
plot(YPred,'Color',[0.6350 0.0780 0.1840],'LineStyle','--','LineWidth',1.5)
legend('YTest','YOutput','Location','northwest')
title(['MAE = ',num2str(MAE_Pred),' , RMSE = ',num2str(RMSE_Pred),' , MAPE = ',num2str(MAPE_Pred),' %'])

here is the test data results :

and here is the prediction data results :

also my data and functions are attached.