function ODEwithEvents
  clc;

# parameters oscillator
  param.m = 1;
  param.b = 0;
  param.c = 1;

# initial values
  x0 = [1;...
        0];

# problem definition
  model = @(t,x) OdeFcn(t,x,param);

# initilaization
  tcum  = [];
  xcum  = [];
  tecum = [];
  xecum = [];
  iecum = [];
  
  tStart = 0;
  tStop  = 10;
  #-----------------------------------------------------------------------------
  # output step size!!!!!!!
  dt   = 1.0e-2;     # ode23/ode45
  #dt   = 1.0001e-2;  # ode15s
  %-----------------------------------------------------------------------------
  tout = (tStart:dt:tStart+dt*(2+fix((tStop-tStart)/dt)))';
  
  tcur    = tStart;
  tbre    = tStop;
  dtevent = 0;

  # calculation
  warning("off","integrate_adaptive:unexpected_termination");
  while tcur < tStop
    
    CurOdeEventFcn = @(t,x) OdeEventFcn(t,x,dtevent);
    opt            = odeset('Events',CurOdeEventFcn,'AbsTol',1e-6,'RelTol',1e-6);
    
    #---------------------------------------------------------------------------
    [t, x, te, xe, ie] = ode23(model,tout,x0,opt);
    #[t, x, te, xe, ie] = ode15s(model,tout,x0,opt);
    #---------------------------------------------------------------------------
    disp(ie)
    disp(te)
    disp('---')
    
    if ~isempty(te)
      tcum = [tcum;t(1:end-1,:)];
      xcum = [xcum;x(1:end-1,:)];
      tecum = [tecum;te(end,:)];
      xecum = [xecum;xe(end,:)];
      iecum = [iecum;ie(end,:)];
      
      tout = (te(end,:):dt:te(end,:)+dt*(2+fix((tStop-te(end,:))/dt)))';
      tcur = te(end,:);
      
      x0 = xe(end,:)';
      #-------------------------------------------------------------------------
      # event count!!!!!!!
      # dependent on solver 1,2,... (ode23/ode45) or
      #                     0,1,... (ode15s/ode15i?)
      if ie == 1 #0 
        x0(2) = -0.5*x0(2);
      elseif ie == 2 #1
        x0(2) = 2*x0(2);
        dtevent = te;
      endif
      #-------------------------------------------------------------------------
    else
      tcum = [tcum;t];
      xcum = [xcum;x];
      tcur = tcum(end);
    endif
  endwhile
  warning("on","integrate_adaptive:unexpected_termination");
  
# visualization
  figure(1);
  subplot(2,1,1);
  plot(tcum, xcum(:,1), '-s');
  grid on;
  subplot(2,1,2);
  plot(tcum, xcum(:,2), '-s');
  grid on;
  
endfunction


# OdeFcn
function dxdt = OdeFcn(t,x,param)
  
  dxdt    = zeros(2,1);
  dxdt(1) = x(2);
  dxdt(2) = -param.b/param.m * x(2) - param.c/param.m * x(1); 
  
endfunction


# OdeEventFcn
function [value,isterminal,direction] = OdeEventFcn(t,x,dt)
  
  # Event setup - state ans time events
  value      = [x(1);t-(1+dt)];
  isterminal = [1;1];
  direction  = [-1;+1];
  
endfunction