I'm trying to set up a simply program that models a billiard ball on a circular, rotating table. I would like to do something like the matlab package file ballode which models a bouncing ball. I'm having trouble locating when the integration reaches a point r=C (polar coordinates). I'd like to have the ball bounce off the edge of the circular table.

What am I doing wrong?

function billiard

%initial conditions
InitRSpeed = 5;
InitPhiSpeed = 0;
InitR = 3;
InitPhi = 4;
y0 = [InitR; InitPhi; InitRSpeed; InitPhiSpeed];

%Rate of rotation of frame
Omega = 1;

%Time
t0 = 0;
tfinal = 20;
Deltat = 0.005*tfinal;
tspan = t0:Deltat:tfinal;

options = odeset('Events',@events,'RelTol',1e-5,'AbsTol',1e-4,'OutputSel',[1 3],'OutputFcn',@odeprint);

%Run DiffEQ
[t,x,te,xe,ie] = ode45(@f,tspan,y0,options);

% Extract the first and third columns of the solution array
Rposn = x(:,1);
Phiposn = x(:,3);
% Plot the result.

polar(Rposn,Phiposn)
title('Polar plot for motion in 2D rotating reference frame')

function dxdt = f(t,x);
Omega=1;
dxdt = [x(2);2*Omega*x(2)*x(4)+x(1)*Omega^2;x(1)*x(4);-2*Omega*x(2)];

function [value,isterminal,direction] = events(t,x)
% Locate when r=boundary and dr/dt positive
value = (x(1)-200); % Detect r = boundary
isterminal = 1; % Stop the integration
direction = 1; % positive direction only