Finding the roots of a function gotten from using ODE45

I'm terrible with MATLAB and I've run into a problem I just can't figure out. Been working on it the past few hours with no luck. Basically, I have a differential equation that I'm using ode45 to solve, and when I get that, I want to use fsolve to find the root of said function. Is there any way to do this, or am I boned? Thanks in advance.

[Star Strider]

If I understand you correctly, you want to find the zero-crossings of your integrated function. That requires inelegant searching (and perhpas sequential searching) using find and then perhaps interpolation to find the most precise value of the root.

Unless you have an analytic expression for your integrated function, you'll not be able to use fsolve.

Answer 1

[Teja Muppirala]

Not as painless as Matt's way, but possible a bit more robust, since it uses the ODE solver's internal zero-finding capabilities (using ODE "Events"):

function findodezeros
options = odeset('Events',@myevents);
vd = @(t,y,mu) [y(2,:);2*(1-y(1,:).^2).*y(2,:)-y(1,:)];
TE = 0;
YE = [2 0];
TEnew = nan;

while ~isempty(TEnew)
    [TOUT,YOUT,TEnew,YEnew]=ode45(vd,[TE(end) 20],YE(end,:),options);
    TE = [TE; TEnew];
    YE = [YE; YEnew];
    plot(TOUT,YOUT); hold all;
end
TE(1) = []
plot(TE,0,'x');
function [value,isterminal,direction] = myevents(t,y)
value = y;
isterminal = ones(size(y));
direction = zeros(size(y));

Answer 2

[Matt Fig]

You can do it fairly painlessly. Here is an example of finding the zeros of the solution to the Van der Pol equation with mu=2:

function [U,V] = ode_examp()
% Solve the Van der Pol DE for mu = 2, then find the zeros.
vd = @(t,y,mu) [y(2,:);2*(1-y(1,:).^2).*y(2,:)-y(1,:)];
[t,y]=ode45(vd,[0 20],[2 0]);
plot(t,y(:,1),'r',t,y(:,2),'b');
f = @(x) interp1(t,y(:,1),x,'spline');
g = @(x) interp1(t,y(:,2),x,'spline');

for ii = 2:2:19  % based on the plot.
    try
        V(ii) = fzero(f,ii);
        U(ii) = fzero(g,ii);
    catch
    end
end

hold on
V = unique(V);
U = unique(U);
L = plot(V,zeros(size(V)),'ob',U,zeros(size(U)),'or');
set(L,{'markerfacecol'},{'r';'b'},'linewidth',2)