function [t,y,y_theta,x] = mtt_stime(system_name,theta,free);
  ## usage: [t,y,y_theta] = mtt_stime(system_name,theta);
  ##
  ## Simulate system with name system_name and parameter vector theta
  ## The order of components in theta is determined in system_numpar.txt:
  ## y_theta contains the corresponding sensitivity functions
  ## Assumes system generated by the sBG approach
  ## Copyright (C) 1999 by Peter J. Gawthrop

  ## 	$Id$	

  ## Simulate using mtt-generated function

  N = length(theta);

  eval(sprintf("[nx,ny,nu,nz,nyz] = %s_def;", system_name));
  if nargin<3
    free = 1;
  endif
  
  y_theta = [];

  if length(free)==0
    free=1;			# Make the loop happen once to get y and X
  endif
  
  [n,m]  = size(free);
  if m==1
    free = free';
  endif
  
  for i=free
    args=sprintf("%i",i);
    for j=1:length(theta)
      args = sprintf("%s %g",args, theta(j));
    endfor

    ## Run system and replace NaN by 1e30 -- easier to handle
    command = sprintf("./%s_ode2odes.out %s | sed \'s/NAN/1e30/g\' >mtt_data.dat\n", \
    system_name, args);
   system(command);

    ## Retrieve data
    load -force mtt_data.dat
    y_theta = [y_theta mtt_data(:,3:2:1+ny)];
  endfor

  ## System data
  [n,m]=size(mtt_data);
  t = mtt_data(:,1);
  y = mtt_data(:,2:2:ny);
  x = mtt_data(:,3+ny:m);

endfunction