if !struct_contains(p_c,"delta_ol")
    p_c.delta_ol = 0.5;	# OL sample interval
  endif
  
  if !struct_contains(p_c,"T")
    p_c.T = 10;			# Last time point.
  endif

  if !struct_contains(p_c,"Iterations")
    p_c.Iterations = 5;		# Number of interations, total =T*Iterations
  endif

  if !struct_contains(p_c,"augment")
    p_c.augment = 0;		# Augment basis funs with constant
  endif
  
  if !struct_contains(p_c,"integrate")
    p_c.integrate = 0;		

    [ol_zeros, ol_poles] = sys2zp(sys)
    cl_poles = eig(A - B*k_x)
    t_max = 1/min(abs(cl_poles))
    t_min = 1/max(abs(cl_poles))
  endif

  ## Initial control U
  U = zeros(p_c.n_U,1);	
  U = zeros(p_c.n_U,1)	


  ## Short sample interval
  dt = p_c.delta_ol/p_c.N;

  ## Observer design
  G = eye(n_x);		# State noise gain 
  sigma_x = eye(n_x);		# State noise variance
  Sigma = p_o.sigma*eye(n_y);	# Measurement noise variance
  Sigma = p_o.sigma*eye(n_y)	# Measurement noise variance
  
  if strcmp(p_o.method, "intermittent")
    Ad = expm(A*p_c.delta_ol);		# Discrete-time transition matrix
    if (ControlType==2)		# 
      [L, M, P, obs_poles] = dlqe(Ad,G,C,sigma_x,Sigma);
    else
      L = zeros(n_x,n_y);

  u = [];
  t = [];
  y_e = [];
  t_e = [];
  e_e = [];
  tick = time;
  i=0;

  for j=1:4
  for j=1:p_c.Iterations
    for k=1:I
      tim=time;			# Timing
      i++;
      if Simulate		# Exact simulation 
	t_sim = [1:p_c.N]*dt;	# Simulation time points
	[yi,ui,xsi] = ppp_ystar(A,B,C,D,x,p_c.A_u,U,t_sim); # Simulate
	x = xsi(:,p_c.N);	# Current state (for next time)

	y = [y;yi'];
	u = [u;ui'];
      else
	t = [t;t_i];
	y = [y;y_i'];
	u = [u;u_i'];
      endif
      

      if strcmp(p_o.method, "intermittent")
	y_e = [y_e; y_new'];
	e_e = [e_e; e_est'];
	t_e = [t_e; t_i];
      endif
      if !Simulate