function event_first_step y0 = [1; 0; 0; 0.3]; eFunc = @(t,y)gstop(t,y,y0); opts = odeset('events',eFunc); tspan=linspace(0, 2*pi, 10); warning ("off", "integrate_adaptive:unexpected_termination", "local"); [t,y,te,ye,ie] = ode45(@(t,y)twobody(t,y,y0),tspan,y0,opts); ie tfinal = te(end) yfinal = ye(end,:) figure;plot(y(:,1),y(:,2),'-',0,0,'ro') axis([-.1 1.05 -.35 .35]) end function ydot = twobody(t,y,y0) r = sqrt(y(1)^2 + y(2)^2); ydot = [y(3); y(4); -y(1)/r^3; -y(2)/r^3]; end function [val,isterm,dir] = gstop(t,y,y0) d = y(1:2)-y0(1:2); v = y(3:4); val = d'*v; isterm = 1; dir = 1; end