/* -*- MaTX -*-
 *
 *  NAME
 *      ltitr() - Time response of discrete-time linear system
 *
 *  SYNOPSIS
 *      X = ltitr(A, B, U)
 *         Array X;
 *         Matrix A,B;
 *         Array U;
 *
 *      X = ltitr(A, B, U, x0)
 *         Array X;
 *         Matrix A,B;
 *         Array U;
 *         Matrix x0;
 *
 *  DESCRIPTION
 *      ltitr() calculates the time response of the discrete-time system:
 *
 *         x[n+1] = Ax[n] + Bu[n]
 *
 *      to input sequence U, where Rows(U) = Rows(u).  Each column of U 
 *      corresponds to a time point.
 *
 *      ltitr() returns the Rows(x)-by-Cols(U) state-history matrix X.
 *
 *      ltitr(..., x0) uses x0 to specify the initial conditions. 
 *
 *  SEE ALSO
 *      ltifr
 */ 

Func Matrix ltitr(A, B, U_, x0, ...)
	Matrix A;
	Matrix B;
	Matrix U_;
	Matrix x0;
{
	Integer i, n, N;
	Matrix X, U, u, x;
	String msg;

	error(nargchk(3, 4, nargs, "ltitr"));

	if (length((msg = abcdchk(A, B))) > 0) {
		error("ltitr(): " + msg);
	}
	
	if (Rows(U_) > Cols(U_)) {
		U = trans(U_);
	} else {
		U = U_;
	}

	n = Cols(A);
	N = length(U);
	
	if (nargs == 3) {
		x0 = Z(n,1);
	}
	
	if (iscomplex(A) || iscomplex(B) || iscomplex(x0)) {
		X = (Z(n,N),:);
	} else {
		X = Z(n,N);
	}

	x = x0;
	for (i = 1; i <= N; i++) {
		X(:,i) = x;
		u = U(:,i);
		x = A*x + B*u;
	}

	if (Rows(U_) > Cols(U_)) {
		return trans(X);
	} else {
		return X;
	}
}