--- a/main/control/devel/lsim2.m	Sun Sep 23 12:05:54 2012 +0000
+++ b/main/control/devel/lsim2.m	Sun Sep 23 13:04:26 2012 +0000
@@ -1,4 +1,4 @@
-## Copyright (C) 2009, 2010, 2011   Lukas F. Reichlin
+## Copyright (C) 2009, 2010, 2011, 2012   Lukas F. Reichlin
 ##
 ## This file is part of LTI Syncope.
 ##
@@ -60,7 +60,7 @@
 
 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
 ## Created: October 2009
-## Version: 0.3
+## Version: 0.4
 
 % function [y_r, t_r, x_r] = lsim (sys, u, t = [], x0 = [], method = "zoh")
 function [y_r, t_r, x_r] = lsim (varargin)
@@ -69,19 +69,24 @@
     print_usage ();
   endif
 
-  sys_idx = cellfun (@isa, varargin, {"lti"});                          # look for LTI models
-  sys_cell = cellfun (@ss, varargin(sys_idx), "uniformoutput", false);  # convert to state-space
+  sys_idx = find (cellfun (@isa, args, {"lti"}));                   # look for LTI models, 'find' needed for plot styles
+  sys_cell = cellfun (@ss, args(sys_idx), "uniformoutput", false);  # convert to state-space
 
   if (! size_equal (sys_cell{:}))
     error ("lsim: models must have equal sizes");
   endif
 
-  vec_idx = find (cellfun (@is_real_matrix, varargin));
-  n_vec = length (vec_idx);
-  n_sys = length (sys_cell);
+  mat_idx = find (cellfun (@is_real_matrix, args));                 # indices of matrix arguments
+  n_mat = length (mat_idx);                                         # number of vector arguments
+  n_sys = length (sys_cell);                                        # number of LTI systems
 
 
+  ## function [y, x_arr] = __linear_simulation__ (sys_dt, u, t, x0)
+  
+  [y, x] = cellfun (@__linear_simulation__, sys_dt_cell, {u}, t, {x0}, "uniformoutput", false);
 
+
+%{
   if (! isa (sys, "ss"))
     sys = ss (sys);                             # sys must be proper
   endif
@@ -148,9 +153,82 @@
     error ("lsim: input vector u must have %d columns", m);
   endif
 
+%}
+
+  if (nargout == 0)                             # plot information
+    [p, m] = size (sys_cell{1});
+    style_idx = find (cellfun (@ischar, varargin));
+    str = "Linear Simulation Results";
+    outname = get (sys_cell{end}, "outname");
+    outname = __labels__ (outname, "y");
+    colororder = get (gca, "colororder");
+    rc = rows (colororder);
+
+    for k = 1 : n_sys                                   # for every system
+      if (k == n_sys)
+        lim = numel (args);
+      else
+        lim = sys_idx(k+1);
+      endif
+      style = args(style_idx(style_idx > sys_idx(k) & style_idx <= lim));
+      if (isempty (style))
+        color = colororder(1+rem (k-1, rc), :);
+        style = {"color", color};   
+      endif
+      if (ct_idx(k))                                    # continuous-time system                                           
+        for i = 1 : p                                   # for every output
+            subplot (p, 1, i);
+            plot (t{k}, y{k}(:, i), style{:});
+            hold on;
+            grid on;
+            if (k == n_sys)
+              axis tight
+              ylim (__axis_margin__ (ylim))
+              ylabel (outname{i});
+              if (i == 1)
+                title (str);
+              endif
+            endif
+          endfor
+        endfor
+      else                                              # discrete-time system
+        for i = 1 : p                                   # for every output
+            subplot (p, 1, i);
+            stairs (t{k}, y{k}(:, i), style{:});
+            hold on;
+            grid on;
+            if (k == n_sys)
+              axis tight;
+              ylim (__axis_margin__ (ylim))
+              ylabel (outname{i});
+              if (i == 1)
+                  title (str);
+              endif
+            endif
+          endfor
+        endfor
+      endif
+    endfor
+    xlabel ("Time [s]");
+    if (p == 1 && m == 1)
+      legend (sysname)
+    endif
+    hold off;
+  endif
+  
+endfunction
+
+
+function [y, x_arr] = __linear_simulation__ (sys_dt, u, t, x0)
+
+  [A, B, C, D] = ssdata (sys_dt);
+  [p, m] = size (D);                            # number of outputs and inputs
+  n = rows (A);                                 # number of states
+  len_t = length (t);
+
   ## preallocate memory
-  y = zeros (trows, p);
-  x_arr = zeros (trows, n);
+  y = zeros (len_t, p);
+  x_arr = zeros (len_t, n);
 
   ## initial conditions
   if (isempty (x0))
@@ -162,46 +240,13 @@
   x = reshape (x0, [], 1);                      # make sure that x is a column vector
 
   ## simulation
-  for k = 1 : trows
+  for k = 1 : len_t
     y(k, :) = C * x  +  D * u(k, :).';
     x_arr(k, :) = x;
     x = A * x  +  B * u(k, :).';
   endfor
 
-  if (nargout == 0)                             # plot information
-    str = ["Linear Simulation Results of ", inputname(1)];
-    outname = get (sys, "outname");
-    outname = __labels__ (outname, "y_");
-    if (discrete)                               # discrete system
-      for k = 1 : p
-        subplot (p, 1, k);
-        stairs (t, y(:, k));
-        grid ("on");
-        if (k == 1)
-          title (str);
-        endif
-        ylabel (sprintf ("Amplitude %s", outname{k}));
-      endfor
-      xlabel ("Time [s]");
-    else                                        # continuous system
-      for k = 1 : p
-        subplot (p, 1, k);
-        plot (t, y(:, k));
-        grid ("on");
-        if (k == 1)
-          title (str);
-        endif
-        ylabel (sprintf ("Amplitude %s", outname{k}));
-      endfor
-      xlabel ("Time [s]");
-    endif
-  else                                          # return values
-    y_r = y;
-    t_r = t;
-    x_r = x_arr;
-  endif
-
 endfunction
 
 
-## TODO: add test cases
\ No newline at end of file
+## TODO: add test cases

--- a/main/control/devel/multiplot2.m	Sun Sep 23 12:05:54 2012 +0000
+++ b/main/control/devel/multiplot2.m	Sun Sep 23 13:04:26 2012 +0000
@@ -94,7 +94,7 @@
 
 % Step Response of Closed Loop
 figure (2)
-step2 (T, Tr4, Tr2)
+step (T, Tr4, Tr2)
 
 %{    
 % Step Response of Closed Loop

--- a/main/control/inst/__time_response__.m	Sun Sep 23 12:05:54 2012 +0000
+++ b/main/control/inst/__time_response__.m	Sun Sep 23 13:04:26 2012 +0000
@@ -166,16 +166,15 @@
         color = colororder(1+rem (k-1, rc), :);
         style = {"color", color};   
       endif
-      discrete = ! ct_idx(k);
-      if (discrete)                                     # discrete-time system
+      if (ct_idx(k))                                    # continuous-time system                                           
         for i = 1 : p                                   # for every output
           for j = 1 : cols                              # for every input (except for initial where cols=1)
             subplot (p, cols, (i-1)*cols+j);
-            stairs (t{k}, y{k}(:, i, j), style{:});
+            plot (t{k}, y{k}(:, i, j), style{:});
             hold on;
             grid on;
             if (k == n_sys)
-              axis tight;
+              axis tight
               ylim (__axis_margin__ (ylim))
               if (j == 1)
                 ylabel (outname{i});
@@ -186,20 +185,15 @@
             endif
           endfor
         endfor
-      else                                              # continuous-time system
+      else                                              # discrete-time system
         for i = 1 : p                                   # for every output
           for j = 1 : cols                              # for every input (except for initial where cols=1)
             subplot (p, cols, (i-1)*cols+j);
-            ##if (n_sys == 1 && isstable (sys_cell{1}))
-            ##  plot (t{k}, y{k}(:, i, j), style{:}, [t{k}(1), t{k}(end)], repmat (yfinal(i,j), 1, 2));
-            ##  ## TODO: plot final value first such that its line doesn't overprint the response
-            ##else
-            plot (t{k}, y{k}(:, i, j), style{:});
-            ##endif
+            stairs (t{k}, y{k}(:, i, j), style{:});
             hold on;
             grid on;
             if (k == n_sys)
-              axis tight
+              axis tight;
               ylim (__axis_margin__ (ylim))
               if (j == 1)
                 ylabel (outname{i});