function [sys,x0,str,ts] = stconcon(t,x,u,flag,tsamp) 
  switch flag,
  case 0,
    [sys,x0,str,ts]=mdlInitializeSizes(tsamp);
  case 2,
    sys=mdlUpdate(t,x,u);
  case 3,
    sys=mdlOutputs(t,x,u);
  case 9,
    sys=mdlTerminate(t,x,u);
  otherwise
    error(['Unhandled flag = ',num2str(flag)]);
end

function [sys,x0,str,ts]=mdlInitializeSizes(tsamp)
  sizes = simsizes;
  sizes.NumContStates  = 0;
  sizes.NumDiscStates  = 0;
  sizes.NumOutputs     = 0;
  sizes.NumInputs      = 5;
  sizes.DirFeedthrough = 1;
  sizes.NumSampleTimes = 1;
  sys = simsizes(sizes);
  x0  = [];
  str = [];
  ts  = [tsamp 0];
  % Controller initialisation
  Denominator = '[.1 .7 .2]';
  Numerator = '[.1 .4 .2]';
  set_param('ad033s/Controller', 'Numerator', Numerator, ...
                                 'Denominator', Denominator);

function sys=mdlOutputs(t,x,u)
  if t>4 % begin adaptive control
    b1 = u(1); b0 = u(2); b = b1*s + b0;
    a2 = u(3); a1 = u(4); a0 = u(5); a = a2*s^2 + a1*s + a0;

    An=[0 0 1 0;0 0 0 1;0 0 0 1;0 0 -a0 -a1];
    Bn=[0;0;0;1];
    Cn=[b0 b1 0 0];
    [BRs,AR]=ss2rmf(An,Bn,eye(4));
    BR=Cn*BRs;

    Q1=1; Q2=1;
    Qn=[Q1 0 0 0;0 Q2 0 0;0 0 0 0;0 0 0 0];
    [Dc,J]=spf(AR'*AR+BRs'*Qn*BRs);
    Df=s^2+1.5*s+3;
    [x0,y0]=axbyc(AR,BR,Dc*Df);
    tt = 0;
    xx = x0 - b * tt;
    yy = y0 + a * tt;

    Denominator = sprintf('[%f %f %f]',xx{2},xx{1},xx{0});
    Numerator = sprintf('[%f %f %f]',yy{2},yy{1}, yy{0});
    set_param('ad033s/Controller', 'Numerator', Numerator, ...
                                   'Denominator', Denominator);
  end
  sys=[];

function sys=mdlUpdate(t,x,u)
  sys = [];

function sys=mdlTerminate(t,x,u)
  sys = [];
