function [t,y,y_theta,x,x_last] = mtt_stime(system_name,theta,free);
  ## usage: [t,y,y_theta,x,x_last] = 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
    file_name = sprintf("%s_input.dat", system_name);
    if exist(file_name)==2	# Then use data from this file ...
      command = sprintf("./%s_ode2odes.out %s < %s | sed \'s/NAN/Inf/g\' >mtt_out_data.dat\n", ...
			system_name, args, file_name);
    else
      command = sprintf("./%s_ode2odes.out %s | sed \'s/NAN/Inf/g\' >mtt_out_data.dat\n", ...
			system_name, args);
    endif
      
   system(command);

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

  ## System out_data
  [n,m]=size(mtt_out_data);
  t = mtt_out_data(:,1);
  y = mtt_out_data(:,2:2:ny);
  x = mtt_out_data(:,3+ny:m);
  x_last = mtt_out_data(n,3+ny:m);
endfunction