function [y_u, Us] = ppp_y_u (A,B,C,D,A_u,u_0,tau)

  ## usage:  y_u = ppp_y_u (A,B,C,D,A_u,u_0,t)
  ##
  ## Computes y_u derivative of y* wrt U
  ## Called by ppp_lin
  ## OBSOLETE

  ###############################################################
  ## Version control history
  ###############################################################
  ## $Id$
  ## $Log$
  ## Revision 1.6  2000/12/27 16:41:05  peterg
  ## *** empty log message ***
  ##
  ## Revision 1.5  1999/05/31 01:58:01  peterg
  ## Now uses ppp_extract
  ##
  ## Revision 1.4  1999/05/12 00:10:35  peterg
  ## Modified for alternative (square) A_u
  ##
  ## Revision 1.3  1999/05/03 23:56:32  peterg
  ## Multivariable version - tested for N uncoupled systems
  ##
  ## Revision 1.2  1999/05/03 00:38:32  peterg
  ## Changed data storage:
  ## y_u saved as row vector, one row for each time, one column for
  ## each U
  ## y_x saved as row vector, one row for each time, one column for
  ## each x
  ## W* saved as row vector, one row for each time, one column for
  ## each element of W*
  ## This is consistent with paper.
  ##
  ## Revision 1.1  1999/04/29 06:02:43  peterg
  ## Initial revision
  ##
  ###############################################################



  ## Check argument dimensions
  [n_x,n_u,n_y] = abcddim(A,B,C,D);
  if (n_x==-1)
    return
  endif

  [n,m] = size(A_u);		# Size of composite A_u matrix
  N = m;			# Number of U* functions per input
  
  y_u = [];			# Initialise
  Us = [];
  
#   for input=1:n_u		# Do for each system input
#     a_u = ppp_extract(A_u,input); # Extract the relecant A_u matrix
#     for i=1:N			# Do for each input function U*_i
#       C_u = zeros(1,N); C_u(i) = 1; # Create C_u for this U*_i
#       b = B(:,input);		# B vector for this input
#       d = D(:,input);		# D vector for this input
#       [y,u] = ppp_transient (t,a_u,C_u,u_0,A,b,C,d); # Compute response for this input
#       y_u = [y_u y'];		# Save y_u (y for input u)  with one row for each t.
#       Us = [Us u'];		# Save u (input)  with one row for each t.
#     endfor
#   endfor
  i_U = 0;
  x_0 = zeros(n_x,1);		# This is for x=0
  for input=1:n_u		# Do for each system input
    a_u = ppp_extract(A_u,input); # Extract the relevant A_u matrix
    for i=1:N			# Do for each input function U*_i
      dU = zeros(N*n_u,1);
      dU(++i_U) = 1;		# Create dU/dU_i 
      [ys,us] = ppp_ystar (A,B,C,D,x_0,a_u,dU,tau); # Find ystar and ustar
      y_u = [y_u ys'];		# Save y_u (y for input u)  with one row for each t.
      Us = [Us us'];		# Save u (input)  with one row for each t.
    endfor
  endfor

endfunction