Passing time step from ODE solver (ode23tb) to ODE function

I am currently utilizing MATLAB's ODE solver ode23tb to advance a system of equations with the structure:

[t, X] = ode23tb(@myfun, [t1 t2], X0)

With the ODE system defined as:

function dxdt = myfun(t, x)
    % Some operations to figure out the time rate of x based on current values of x, stored ultimately in a variable 'blah'
    dxdt = blah;
end

Here is the catch. The evaluation of 'blah' depends on all previous steps taken. (My governing equation involves fractional derivative, which fundamentally has a convolution integral ... without getting into too much detail.) Therefore I need to find an array of the successful steps taken by ode23tb.

So far, the only way I have been able to make this work is to create a modified version of ode23tb and globalize tout. (I also had to place a copy of the private function folder in the local directory.) By checking tout inside myfun, I can get a list of the steps that ode23tb stored as time history of the ODE advancement.

However, declaring global variable is obviously a bad practice. I also need to do parameter sweep eventually and run myfun multiple times with different input parameters. Yet, the global declaration creates a mess that can only be cleaned with the horrible clear all command -- which makes it challenging to run parameter sweep with a loop structure.

What would be the proper way to retrieve the tout information from ode23tb without generating a monster comprised of global and clear all?

Thanks in advance for the useful feedbacks!

[Zhiren]

Right after posting the question, I came up with a weird method: saving the variable 'tout' in an intermediate .mat file whenever ode23tb gets new step and then loading it every time myfun needs 'tout'.

Pretty sure this is not the best way to do this. Looking forward to further improving my approach!

Answer 1

[Jan]

You can use the OutputFcn to store the successful steps in a persitent variable, which can be forwarded to the function to bbe integrated. This is ways faster than writing to disk.

function status = myOutputFcn(t,y,flag)
persistent T_local, T_count

switch char(flag)
    case ''
        T_count = T_count + 1;
        if T_count > numel(T_local)
            T_local(T_count + 1000) = 0;
        end
        T_local(t_count) = t;
        status = 0;
        
    case {'init', 'done'}
        T_local = zeros(1, 1000);
        T_count = 0;
        status  = 0;
        
    case 'get'
        status = T_local(1:T_count);
end
end

Now the funciton to be integrated cann call:

formerT = myOutputFcn([], [], 'get')

[Zhiren]

This is clever. Thanks!

My saving to disk approach indeed slowed down the process quite significantly. Glad to find a faster alternative -- and so promptly!

Answer 2

[Steven Lord]

The evaluation of 'blah' depends on all previous steps taken.

That sounds like you don't have an ordinary differential equation but a delay differential equation. In that case look at dde23.

[Zhiren]

Oh interesting. I am actually unfamiliar with the concept of DDE. Will certainly look into that. Thanks for the pointer!