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