If you use ode45 then a chunk of your code becomes moot - you just need the initial conditions y0 and time span [t0 tfinal]. Your equations are encoded into a function that you pass to the solver. Note the use of anonymous functions here to pass the extra parameters to the ode and event functions.
%Inputs
tspan = [50 55];
y0 = [-A/(M*(2*pi*F)^2) 0 0 0];
opts = odeset('Events',@(t,y) impactEvents(t,y,d,w,r));
[t,y,te,ye,ie] = ode23(@(t,y) odefun(t,y,A,F,M,d,w,r,m),tspan,y0,opts);
function yxzv = odefun(t,y,A,F,M,d,w,r,m)
dydt = y(2);
dxdt = (A*cos(2*pi*F*t))/M;
dzdt = y(4);
dvdt = 0;
yxzv = [dydt; dxdt; dzdt; dvdt];
end
The impact is controlled by the event function, which tells you when a certain expression is equal to zero. Due to the way your impact condition is expressed (one quantity larger than another), it needs to be rearranged to fit that format (one quantity equal to zero). Here is a quick first pass at what the event function might look like:
function [value,isterminal,direction] = impactEvents(t,y,d,w,r)
value = y(1)-y(3)-d-2*w+2*r; %value that describes the event
isterminal = 1; %should the integration stop when an event occurs?
direction = 1; %From what direction will the value approach zero, from above or below?
end
One snag I noticed immediately is that the way the impact condition is expressed above (which might be wrong), no events ever occur. So you might want to double check the condition, equations, and/or direction.
>> any((y(:,1)-y(:,3))>(d+2*w-2*r))
ans =
logical
0
Once you get that working, you'll want to add some more code to restart the integration after an event occurs. The velocity after the impact will be the new initial condition for the velocity, and then the integration proceeds until it finishes or another event occurs. For an example of how to stop/restart based on events type edit ballode - this is an example of a bouncing ball.
