/* -*- MaTX -*-
 *
 *  NAME
 *      dstep() - Step response of discrete-time linear systems
 *
 *  SYNOPSIS
 *      {Y,X} = dstep(A,B,C,D,iu,N)
 *         Matrix Y,X;
 *         Matrix A,B,C,D;
 *         Integer iu,N;
 *
 *  DESCRIPTION
 *       dstep(A,B,C,D,iu,N) calculates the response of the system:
 *
 *      	x[n+1] = Ax[n] + Bu[n]
 *      	y[n]   = Cx[n] + Du[n]
 *
 *       to an unit step applied to the iu'th input.
 *
 *       dimpulse() returns the Rows(y)-by-N output-history matrix Y and
 *       the state-history matrix X.
 *
 *       If iu = 0, dstep() returns matrices
 *
 *       [[Y for 1st input]          [[X for 1st input]
 *        [Y for 2nd input]           [X for 2nd input]
 *        [...............]           [...............]
 *        [Y for m'th input]]  and    [Y for m'th input]].
 *
 *  SEE ALSO
 *      step and dimpuse
 */

Func List dstep(A, B, C, D, iu, N)
	Matrix A, B, C, D;
	Integer iu, N;
{
	Integer i,m;
	Matrix X,Y,XX,YY;
	String msg;

	if (length((msg = abcdchk(A,B,C,D))) > 0) {
		error("dstep(): " + msg);
	}

	if (iu == 0) {
		m = Cols(B);
		YY = Z(Rows(C)*m,N);
		XX = Z(Rows(A)*m,N);
		if (version() == 4) {
			for (i = 1; i <= m; i++) {
				{Y,X} = dlsim(A, B(:,i), C, D(:,i), ONE(1,N));
				YY(i-1,0,Y) = Y;
				XX(i-1,0,X) = X;
			}
		} else {
			for (i = 1; i <= m; i++) {
				{Y,X} = dlsim(A, B(:,i), C, D(:,i), ONE(1,N));
				YY(i,1,Y) = Y;
				XX(i,1,X) = X;
			}
		}
	} else {
		{YY,XX} = dlsim(A, B(:,iu), C, D(:,iu), ONE(1,N));
	}
	return {YY,XX};
}