How to plot the model prediction error as a function of forecast horizon?

I would like to be able to identify a model from some data, say the first 1000 time samples, and then use the remaining data (800 points) to evaluate the model in the following way. Fit the model initial condition to give the best response for Np points, and then evaluate the model for the next Nf points, where Np are past points, and Nf are future points. Then determine the error between the model and the actual data for the Nf future points. Repeat by shifting the starting point, till all available data is used. The plot error distribution for values of Nf, ie, Nf = 5, 15, 30, 60.

Following an example for n4sid

close all; clear all;
N = 1800; K = 0.5; Knl=0.2;
rng('default');
w = randn(N,1);
z = zeros(N,1); u = zeros(N,1); y = u;
e = randn(N,1);
v = filter([1 0.5],[1 1.5 0.7],e);
z(1) = 0; z(2) = 0; u(1) = 0; u(2) = 0; y(1) = 0; y(2) = 0;
for k = 3:N
    u(k-1) = -K*y(k-2)+w(k)-Knl*y(k-2)^2;
    u(k-1) = -K*y(k-1)+w(k);
    z(k) = 1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);
    y(k) = z(k) + .8*v(k);
end
dat = iddata(y, u, 1);

optCVA = n4sidOptions('N4weight','CVA');
optSSARX = n4sidOptions('N4weight','SSARX');

sysCVA = n4sid(dat, 2, optCVA);
sysSSARX = n4sid(dat, 2, optSSARX);

compare(dat, sysCVA, sysSSARX);

But, how do I calculate the prediction starting from a number of points, and then calculate the error, and then display the error distribution for different forecast horizons?