new file mode 100644
--- /dev/null
+++ b/scripts/control/base/lsim.m
@@ -0,0 +1,99 @@
+## Copyright (C) 1996 Auburn University.  All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+## for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} lsim (@var{sys}, @var{u}, @var{t}@{,@var{x0}@})
+## Produce output for a linear simulation of a system
+##
+## Produces a plot for the output of the system, sys.
+##
+## U is an array that contains the system's inputs.  Each column in u
+## corresponds to a different time step.  Each row in u corresponds to a
+## different input.  T is an array that contains the time index of the
+## system.  T should be regularly spaced.  If initial conditions are required
+## on the system, the x0 vector should be added to the argument list.
+##
+## When the lsim function is invoked with output parameters:
+## [y,x] = lsim(sys,u,t,[x0])
+## a plot is not displayed, however, the data is returned in y = system output
+## and x = system states.
+## @end deftypefn
+
+## Author: David Clem
+## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu>
+## Created: July 1995
+## modified by John Ingram for system format August 1996
+
+function [y, x] = lsim (sys, u, t, x0)
+
+  if((nargin < 3)||(nargin > 4))
+    usage("[y,x] = lsim(sys,u,t[,x0])");
+  endif
+
+  if(!is_struct(sys))
+    error("sys must be in system data structure");
+  endif
+
+  sys = sysupdate(sys,"ss");
+
+  [ncstates, ndstates, nin, nout] = sysdimensions(sys);
+  [a,b,c,d] = sys2ss(sys);
+
+  if (nargin == 3)     x0 = zeros(columns(a),1);        endif
+
+  if(rows(u) ~= length(t))
+    error("lsim: There should be an input value (row) for each time instant");
+  endif
+  if(columns(u) ~= columns(d))
+    error("lsim: U and d should have the same number of inputs");
+  endif
+  if(columns(x0) > 1)
+    error("lsim: Initial condition vector should have only one column");
+  endif
+  if(rows(x0) > rows(a))
+    error("lsim: Initial condition vector is too large");
+  endif
+
+  Ts = 0;
+  t(2)-t(1);
+  u=u';
+  n = max(size(t));
+  for ii = 1:(n-1)
+
+    ## check if step size changed
+    if (t(ii+1) - t(ii) != Ts)
+      Ts = t(ii+1) - t(ii);
+      ## [F,G] = c2d(a,b,Ts);
+      dsys = c2d(sys, Ts);
+      [F,G] = sys2ss(dsys);
+    endif
+
+    x(:,ii) = x0;
+    x0 = F*x0 + G*u(:,ii);
+  endfor
+
+  ## pick up last point
+  x(:,n) = x0;
+
+  y = c*x + d*u;
+  if(nargout == 0)
+   plot(t,y);
+   y=[];
+   x=[];
+  endif
+endfunction