How to store all solutions of an ODE on Matlab for multiple values of a parameter

Question

I would like to solve an ODE for multiple values of the parameter p and most importantly, save all the solutions for all the different values.

Till now, I have been using this:

p = -200:+1:300;
time = 0:.01:10;
y0 = [0 0 0 0 0 0 0 0]; 
y = NaN(length(time),length(y0),length(p));

for i=1:length(p)
    [t,y(:,:,i)] = ode45(@myode,time,y0,[],p(i)); 
end

but it has the t predefined, which is not supposed to be.

What I see as a problem, is that I cannot store all y values for all times t and all values of p in a matrix because I cannot use the variable t before the loop. If i use the variable time instead, I will not be able to take advantage of the ODE45 integration which uses its own dt intervals, dependent the nonlinearities it will encounter.

A potential solution I could think of is the following:

p = -200:+1:300;
time = [0 10];
y0 = [0 0 0 0 0 0 0 0];
y = NaN(1,length(y0),length(p));

for i=1:length(p)
    [t,x] = ode45(@ode,time,y0,[],p(i)); 

    y(1:length(t),:,i) = x(:,:);
end

In this case, it manages to store all the values, but, the steps of integration changes from one loop to the other. As a result, for the cases with shortest length(t), since they all coexist in the same matrix, the matrix places that rest, are 0. I understand that I probably need to use interpolation to find the values that have been replaced automatically by 0.

I am using the following line within the loop, but it does not work:

y(1:length(t),:,i) = interp1(x,time)

since I want to interpolate using the tspan, but without success. I am guessing I am doing something wrong with the interpolation.

I would really appreciate any suggestion!!

Answer 1

I think you have a problem due to a misunderstanding of how ODE45 works. If I'm understanding your problem correctly, you want to plot the behavior at certain times, but you are concerned that you can't do this with ODE45. This is a misunderstanding of how ODE45 works. ODE45 uses adaptive time-stepping, and will never use the time array you give it. It will use its own adaptive time-stepping procedure with a local error control mechanism and then use a high order interpolation function to interpolate to the time slices that you desire. Therefore, you can choose the time slices ahead of time regardless of the error control mechanism, as the high order interpolation is designed to have similar error behavior as the DPRK45 method that matlab uses in ODE45. You can refer to their documentation for more information.