How do I forecast an ARIMA model?

I want to forecast a lot of series of 4-6 datapoints with ARIMA (I know this isn't much, but I don't have more data). For each serie I use the code below.

A few series:

• 0 0 0 1 0
• 0 0 0 1 1
• 0 2 2 2 2
• 0 NaN 16 0 16
• 0 8 0 3 3
• 0 4 0 1 1
• 0 0 1 1 1

I only get a forecast (Y=8) for [0 NaN 16 0 16]. For the rest Y=NaN. When I delete the Lbqtest (and only perform the else), I get forecasts for all series. But strange numbers for:

• 0 0 0 1 1 (Y = -9,00719e+15)
• 0 0 1 1 1 (Y = -4,49779e-08)

This is for I = 1 (so there is non-stationarity). Can someone help me with the autocorrelation (Lbqtest) and/or the strange numbers when I don't perform the Lbqtest (because I don't know how it works with the Lbqtest)?

 function [ Y, D, ARp, MAq, YMSE ] = ARIMAToolbox( Serie )

 Input = Serie;
 IndexSerieNaN = find(isnan(Input) == 1);
 Input(IndexSerieNaN) = [];  %Delete NaN

MaxI = size(Input,2)-2;

%Checking for I(n) process
VarianceMatrix = zeros(MaxI,1);
for I = 1:MaxI
    VarianceMatrix(I,1) = var(diff(Input,I));
end
IMatrix = find(ismember(VarianceMatrix, min(VarianceMatrix))==1);

if size(IMatrix,1) > 1 
    IMatrix = IMatrix(1);
end

if var(Input) <= min(VarianceMatrix)   
    I = 0;
else
    I = IMatrix;
end

%Making Data Stationary
if I > 0
    StationaryInput = diff(Input,I);
else
    StationaryInput = Input;
end

[row1,~] = size(StationaryInput);

WNLagPeriod = size(StationaryInput,2)-1; 
MaxAR = size(StationaryInput,2)-1;
MaxMA = size(StationaryInput,2)-1;
numPeriods = 1; 

%Testing for White-Noise
%Ljung-Box Q-test for residual autocorrelation. H (hypothesis) = LBTestH.
%H=0: no autocorrelation (white noise)
%H=1: rejection of null hypothesis (autocorrelation)
[LBTestH,~,~,~] = lbqtest(StationaryInput, WNLagPeriod);

D = I;
%Start Iterations
if LBTestH == 0
    disp(['The process is white noise at I = ', num2str(I)]);
    Y = NaN;
    ARp = NaN;
    MAq = NaN;
    YMSE = NaN;
else
    %Selecting AR & MA process based on AIC
    Z = zeros(MaxAR + 1, MaxMA + 1);
    TotLogL = zeros(MaxAR + 1, MaxMA + 1);

      for MA = 0:MaxMA
          for AR = 0:MaxAR 
              if AR+MA < size(StationaryInput,2)
                  i = AR + 1;
                  j = MA + 1;

                  Mdl = arima(AR,D,MA); 
                  try %try-catch: If no error, go on
                      [EstMdl,EstParamCov,logL,~] = estimate(Mdl,StationaryInput'); 
                      numParams = sum(any(EstParamCov)); 
                      [AIC,BIC] = aicbic(logL,numParams,row1); %Akaike and Bayesian information criteria
                      Z(i,j) = AIC;
                      TotLogL(i,j) = logL;
                  catch %if error
                      Z(i,j) = NaN;
                      TotLogL(i,j) = NaN;
                  end
              end
          end
      end

      [row,column] = find(ismember(Z, min(min(Z), [], 2))); 
      ARp = row - 1;
      MAq = column - 1;
      Mdl = arima(ARp,D,MAq);
      [FinalMdl,FinalEstParamCov,~,info] = estimate(Mdl,StationaryInput');

      [Y,YMSE] = forecast(FinalMdl,numPeriods);
  end

end