Only single-input, single-output models are handled by rbj.
Use rpem for the multiple-input case.

The estimated parameters are returned in the matrix thm.
The kth row of thm contains
the parameters associated with time k; that is,
they are based on the data in the rows up to and including row k in z.
Each row of thm contains the estimated parameters
in the following order.

thm(k,:) = [b1,...,bnb,c1,...,cnc,d1,...,dnd,f1,...,fnf]

yhat is the predicted value of the output,
according to the current model; that is, row k of yhat contains
the predicted value of y(k) based on all past data.

The actual algorithm is selected with the two arguments adm and adg.
These are described under rarx.

The input argument th0 contains the initial
value of the parameters, a row vector consistent with the rows of thm.
The default value of th0 is all zeros.

The arguments P0 and P are
the initial and final values, respectively, of the scaled covariance
matrix of the parameters. See rarx.
The default value of P0 is 104 times
the unit matrix. The arguments phi0, psi0, phi,
and psi contain initial and final values of the
data vector and the gradient vector, respectively. The sizes of these
depend on the chosen model orders. The normal choice of phi0 and psi0 is
to use the outputs from a previous call to rbj with
the same model orders. (This call could be a dummy call with default
input arguments.) The default values of phi0 and psi0 are
all zeros.

Note that the function requires that the delay nk be
larger than 0. If you want nk = 0,
shift the input sequence appropriately and use nk = 1.