--- a/extra/control-devel/devel/test_slident.m	Mon Mar 19 15:35:53 2012 +0000
+++ b/extra/control-devel/devel/test_slident.m	Mon Mar 19 16:16:49 2012 +0000
@@ -2031,7 +2031,7 @@
 rcond = 0.0;
 tol = -1.0;
 
-[a, b, c, d, q, ry, s, k] = slident (Y, U, nobr, meth, alg, jobd, batch, conct, ctrl, rcond, tol)
+[a, b, c, d, q, ry, s, k, x0] = slident (Y, U, nobr, meth, alg, jobd, batch, conct, ctrl, rcond, tol)
 
 
 ae = [  0.8924   0.3887   0.1285   0.1716
@@ -2077,22 +2077,26 @@
 
 
 
-figure (1)
-plot (Y)
+% figure (1)
+% plot (Y)
 
 P = ss (a, b, c, d, 1);
 
-figure (2)
-lsim (P, U)
-%lsim (P, U, [], x0)  % initial values dependent on realization, IB01BD != IB01CD
+% figure (2)
+% lsim (P, U)
+% lsim (P, U, [], x0)  % initial values dependent on realization, IB01BD != IB01CD
 
 
-[y, t] = lsim (P, U);
+% [y, t] = lsim (P, U);
+[y, t] = lsim (P, U, [], x0);
+
 
 figure (3)
 plot (t, Y, 'b', t, y, 'r')
 legend ('y measured', 'y simulated', 'location', 'southeast')
-axis tight
+%axis tight
+
+% print -depsc2 ident.eps
 
 %n
 %sv

--- a/extra/control-devel/src/slident.cc	Mon Mar 19 15:35:53 2012 +0000
+++ b/extra/control-devel/src/slident.cc	Mon Mar 19 16:16:49 2012 +0000
@@ -561,7 +561,6 @@
 //       ldwork = max (ldwork_a, ldwork_b, ldwork_c)
 
 
-/*
         // arguments in
         char jobx0 = 'X';
         char comuse = 'U';
@@ -571,56 +570,18 @@
         int ldv = max (1, n);
         
         ColumnVector x0 (n);
-        Matrix (ldv, n);
+        Matrix v (ldv, n);
         
         // workspace
         int liwork_c = n;     // if  JOBX0 = 'X'  and  COMUSE <> 'C'
         int ldwork_c;
-        
-        ldwork_c = ldw1 + n*( n + m + l ) + max (5*n, ldw1, min (ldw2, ldw3))
-        
-        
-        // I can't find a definition for parameter 't' 
-        
-C             The length of the array DWORK.
-C             LDWORK >= 2,  if  JOBX0 = 'N'  and  COMUSE <> 'C',  or
-C                           if  max( N, M ) = 0.
-C             Otherwise,
-C             LDWORK >= LDW1 + N*( N + M + L ) +
-C                              max( 5*N, LDW1, min( LDW2, LDW3 ) ),
-C             where, if  COMUSE = 'C',  then
-C             LDW1 = 2,          if  M = 0  or   JOB = 'B',
-C             LDW1 = 3,          if  M > 0  and  JOB = 'D',
-C             LDWa = t*L*(r + 1) + max( N + max( d, f ), 6*r ),
-C             LDW2 = LDWa,       if  M = 0  or  JOB = 'B',
-C             LDW2 = max( LDWa, t*L*(r + 1) + 2*M*M + 6*M ),
-C                                if  M > 0  and JOB = 'D',
-C             LDWb = (b + r)*(r + 1) +
-C                     max( q*(r + 1) + N*N*M + c + max( d, f ), 6*r ),
-C             LDW3 = LDWb,       if  M = 0  or  JOB = 'B',
-C             LDW3 = max( LDWb, (b + r)*(r + 1) + 2*M*M + 6*M ),
-C                                if  M > 0  and JOB = 'D',
-C                r = N*M + a,
-C                a = 0,                  if  JOBX0 = 'N',
-C                a = N,                  if  JOBX0 = 'X';
-C                b = 0,                  if  JOB   = 'B',
-C                b = L*M,                if  JOB   = 'D';
-C                c = 0,                  if  JOBX0 = 'N',
-C                c = L*N,                if  JOBX0 = 'X';
-C                d = 0,                  if  JOBX0 = 'N',
-C                d = 2*N*N + N,          if  JOBX0 = 'X';
-C                f = 2*r,                if  JOB   = 'B'   or  M = 0,
-C                f = M + max( 2*r, M ),  if  JOB   = 'D'  and  M > 0;
-C                q = b + r*L;
-C             and, if  JOBX0 = 'X'  and  COMUSE <> 'C',  then
-C             LDW1 = 2,
-C             LDW2 = t*L*(N + 1) + 2*N + max( 2*N*N, 4*N ),
-C             LDW3 = N*(N + 1) + 2*N + max( q*(N + 1) + 2*N*N + L*N,
-C                                           4*N ),
-C                q = N*L.
-C             For good performance,  LDWORK  should be larger.
+        int t = nsmp;
+   
+        int ldw1_c = 2;
+        int ldw2_c = t*l*(n + 1) + 2*n + max (2*n*n, 4*n);
+        int ldw3_c = n*(n + 1) + 2*n + max (n*l*(n + 1) + 2*n*n + l*n, 4*n);
 
-        
+        ldwork_c = ldw1_c + n*( n + m + l ) + max (5*n, ldw1_c, min (ldw2_c, ldw3_c));
         
         OCTAVE_LOCAL_BUFFER (int, iwork_c, liwork_c);
         OCTAVE_LOCAL_BUFFER (double, dwork_c, ldwork_c);
@@ -632,7 +593,7 @@
 
         // SLICOT routine IB01CD
         F77_XFCN (ib01cd, IB01CD,
-                 (jobx0, comuse, job***,
+                 (jobx0, comuse, jobbd,
                   n, m, l,
                   nsmp,
                   a.fortran_vec (), lda,
@@ -676,7 +637,7 @@
 
         error_msg ("ident", info_c, 2, err_msg_c);
         warning_msg ("ident", iwarn_c, 6, warn_msg_c);
-*/        
+      
         
         // return values
         retval(0) = a;
@@ -689,7 +650,7 @@
         retval(6) = s;
         retval(7) = k;
         
-        //retval(8) = x0;
+        retval(8) = x0;
         //retval(0) = octave_value (n);
         //retval(1) = r;
         //retval(2) = sv;