Create a second output for ‘u’ (ode45 will ignore it during the integration) and then return it using a for loop after the integration finishes.
Try this —
tspan = 0:0.001:15;
y0 = [0.2,0.3];
[t,y] = ode45(@(t,y) odefcn(t,y), tspan, y0);
for k = 1:numel(t)
[dx,u(k)] = odefcn(t(k),y(k,:).');
end
figure
plot(t,y(:,1), 'DisplayName','y_1(t)')
hold on
plot(t,y(:,2), 'DisplayName','y_2(t)')
plot(t,u, 'DisplayName','u(t)')
hold off
grid
legend('Location','best')
% set(gca,'XScale','log')
function [dx,u] = odefcn(t,x)
dx = zeros(2,1);
u=[2 -2]*x(1:2)+1;
dx(1)=x(1)-2*x(2);
dx(2)=-x(2)+u;
end
I changed ‘tspan’ because with a long vector, the details of the initial transient disappear from the plot.
.
