You can select the phase plane output, or you can simply substitute the solved values for ‘t’ and ‘v’ (in my code here) back into the original ODE function to get the derivatives.
This code calls the 'odephas2' output function:
% MAPPING: s = v(1), i = v(2)
ode = @(t, v, beta, mu, gamma) [-beta.*v(1).*v(2) + mu - mu.*v(1); beta.*v(1).*v(2) - (gamma + mu).*v(2)];
beta = rand; % Choose The Correct Values
mu = rand;
gamma = rand;
v_init = rand(2,1);
tspan = linspace(1, 20, 50);
opts = odeset( 'OutputFcn','odephas2'); % Set To Plot Phase Plane
[t, v] = ode45(@(t,v) ode(t, v, beta, mu, gamma), tspan, v_init, opts);
figure(2)
plot(v(:,1), v(:,2)) % Plot The Phase Plane Manually
grid
for k1 = 1:length(t)
derivs(:,k1) = ode(t(k1), v(k1,:), beta, mu, gamma); % Calculate The Derivatives
end
figure(3)
plot(derivs(1,:), derivs(2,:)) % Plot Derivatives
grid
