function [Y,X] = ppp_sm2sr(A,B,C,D,T,u0,x0);
  ## Usage [Y,X] = ppp_sm2sr(A,B,C,D,T,u0,x0);
  ## Computes a step response
  ## A,B,C,D- state matrices
  ## T vector of time points
  ## u0 input gain vector: u = u0*unit step.

  ## Copyright (C) 1999 by Peter J. Gawthrop
  ## 	$Id$	

  [Ny,Nu] = size(D);
  [Ny,Nx] = size(C);

  if nargin<6
    u0 = zeros(Nu,1);
    u0(1) = 1;
  end;

  if nargin<7
    x0 = zeros(Nx,1);
  end;

  [N,M] = size(T);
  if M>N
    T = T';
    N = M;
  end;



  one = eye(Nx);

  Y = zeros(Ny,N);
  X = zeros(Nx,N);

  dt = T(2)-T(1);		# Assumes fixed interval
  expAdt = expm(A*dt);		# Compute matrix exponential
  i = 0;
  expAt = one;

  DoingStep = max(abs(u0))>0;	# Is there a step component?
  DoingInitial = max(abs(x0))>0; # Is there an initial component?
  for t = T'
    i=i+1;
    if Nx>0
      x = zeros(Nx,1);
      if DoingStep
	x = x + ( A\(expAt-one) )*B*u0;
      endif
      if DoingInitial
	x = x + expAt*x0;
      endif
      
      expAt = expAt*expAdt;

      X(:,i) = x;
      if Ny>0
	y = C*x + D*u0;
	Y(:,i) = y;
      endif
    elseif Ny>0
      y = D*u0;
      Y(:,i) = y;
    endif
  endfor


endfunction