matlab保存ode45计算中间值并绘图

假如要用ode45数值积分求解以下方程：
$$y^{{}^{\prime \prime }}=\frac{A}{B}ty$$
将其降阶得到
$$\begin{gathered} y_1^{{}^{\prime }}=y_2 \\ {y}_2^{{}^{\prime }}=\frac{A}{B}t{y}_1 \end{gathered}$$
然后我们希望保存计算的中间值$y_1^{{}^{\prime }}$并绘图，可以这样做：

odefunc.m

function [dydt] = odefunc(t,y,A,B)
%% 要积分的微分方程
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = (A/B)*t.*y(1);

%% 从工作区读取record_dy1数组
record_dy1=evalin('base', 'record_dy1');
record_dy1=[record_dy1 y(2)];
%保存record数组到工作区
assignin('base','record_dy1',record_dy1);

%% 从工作区读取积分时间数组
time=evalin('base', 'time');
%一定要round，消除变步长重复值
time=[time round(t,4)];
assignin('base','time',time);
end

ode45_extend.m

%% 微分方程
A = 1;
B = 2;
tspan = [0 5];
%被积函数dy_1,dy_2初值
y0 = [0 0.01];
%% record_dy1是用来保存要记录的中间值
record_dy1=[];
%time是对应的时间点
time=[];
[t,y] = ode45(@(t,y) odefunc(t,y,A,B), tspan, y0);
%plot(t,y(:,1),'-o',t,y(:,2),'-.')
%% 去除变步长带来的重复值
[B, I] = unique(time, 'first');
%重复值所在索引
duplicated_location=setdiff(1:numel(time), I);
%删除重复值
time(duplicated_location)=[];
record_dy1(duplicated_location)=[];
%% 绘图
plot(time,record_dy1);

运行ode45_extend.m得到结果：