function [theta,Theta,Error,Y] = mtt_identify (system_name,y_s,theta_0,criterion,max_iterations)

  ## usage:  [theta,Theta,Error] = mtt_identify (system_name,theta_0,criterion,max_iterations)
  ##      y_s   actual system output
  ##      theta_0   initial parameter estimate
  ##      criterion convergence criterion
  ##      max_iterations limit to number of iterations

  if nargin<4
    criterion = 1e-5;
  endif
  
  if nargin <5
    max_iterations = 20;
  endif

  alpha = 1.0;
  e_last = 1e20;
  error=1e10;
  theta = theta_0;
  Theta = [];
  Error = [];
  Y = [];
  iterations = 0;
  while abs(e_last-error)>criterion
    theta
    iterations = iterations + 1;
    e_last = error;
    [t,y,y_theta] = mtt_ssimulate(system_name,theta); # Simulate system
    Theta = [Theta theta];	# Save parameters
    Y = [Y y];			# Save output
    E = (y - y_s);		# Error(t)
    error = (E'*E);			# Sum the error
    Error = [Error error];
    theta = theta - alpha*((y_theta'*y_theta)\(y_theta'*E));
  endwhile

endfunction