Using persistent variable to pass out extra parameters using ode45 from the MATLAB ODE suite?

In How do I pass out extra parameters using ODE23 or ODE45 from the MATLAB ODE suite? the "MathWorks Support Team" suggests using a persistent variable to pass out extra parameters using ode45 from the MATLAB ODE suite.

However, what if an integration failed and solver will call your function for the previous time step? The ode45 solver uses an adaptive time step, and sometimes integration can fail, so the solver automatically decreases time step, and "has to go back."

Is "MathWorks Support Team" wrong in their suggestion of using a persistent variable? I see that the 'OutputFcn' method can be used to pass out extra parameters, but I do not know how to use it. Could you give me an outline or example code using the 'OutputFcn' method to detect failed steps / flag etc. to give the correct extra parameters? If I use the 'OutputFcn' method, is a global variable needed?

Solution

Here is an example, slightly modified from the documentation.

Define the ODE function. param will be the variable to be exchanged.

function dy = rigid(t, y)

global param

dy = zeros(3,1);    % a column vector
dy(1) = param + y(2) * y(3);
dy(2) = -y(1) * y(3);
dy(3) = -0.51 * y(1) * y(2);

param = dy(1) * dy(3);

The output function must declare the variable to be exchanged as global. myout will be called after each successful step.

function status = myout(t, y, flag)

global param

% Don't ever halt the integration
status = 0;

if strcmp(flag, 'init')
    % Initialization (apparently useless).
    param = -2;

elseif isempty(flag)
     % Main code to update param.
     param = param * mean(t);

elseif strcmp(flag, 'done')
     % Nothing to do.

end

Call ode45 in this way.

% True initialization.
global param
param = -2;

odeOpts = odeset(  ...
    'RelTol', 1e-6,  ...
    'AbsTol', 1e-7,  ...
    'OutputFcn', @myout);

[time, y] = ode45(@rigid, [0 10], [0 1 1], odeOpts);

plot(  ...
    time, y(:,1), '-',   ...
    time, y(:,2), '-.',  ...
    time, y(:,3), '.');