The issue I have is having to compute the derivative (in real time) of the solution produced by ode45 within the events function.
Some pseudo-code to explain what I'm mean is,
function dx = myfunc(t,x,xevent,ie)
persistent xs
persistent dx2
% xevent is the solution at the event
if isempty(dx2) == 1
dx2 = 0;
end
k = sign(dx2*x(1));
if ie(end) == 1
xs = xevent
elseif ie(end) == 2
xs = xevent
end
dx(1) = x(2);
dx(2) = function of k,x(1),x(2), and xs;
dx2 = dx(2);
end
and
function [value,isterminal,direction] = myeventfcn(t,x)
if dx(2)*x(1)>=0
position = [dx(2)*x(1); dx(2)*x(1)];
isterminal = [1; 1];
direction = [1 ; -1 ];
elseif dx(2)x(1)<0
position = [dx(2)*x(1); dx(2)*x(1)];
isterminal = [1; 1];
direction = [-1 ; 1 ];
end
I know that if I didn't need to use the solution at the event within myfunc I could just compute dx=myfunc(t,x) within my event function and use dx(2), yet since xevent is used within myfunc I can't input xevent.
I know there is a way of inputting constant parameters within the event function, but since it's the solution at the event location, which also changes, I'm not sure how to go about doing that.
My work around to this is to approximate dx(2) using the solution x. What I wanted to know is if it will be satisfactory to just use a finite difference approximation here, using a small fixed step size relative to the step size od45 takes, which is a variable step size.
EDIT: I went ahead and updated my code to reflect what myfunc and myeventfcn would look like. As a note, the reason I have myeventfcn split by the if statement is to know what the direction the event is crossed, since it will update within myfunc.
Also, I need to use dx(2) of the previous successful time step and so that's why I have dx2 defined as a persistent variable. I believe having k=sign(dx(2)*x(1)) is okay to do, since my event depends on dx(2)*x(1), and so I won't be introducing any new discontinuities with the sign function.