a = 0.1; b = 0.004; c = 0.2; d = 0.001; y0 = [100, 30];
Answer.
To the code lines written above, we only have to add
rhs = @(t, y) [a*y(1) - b*y(1)*y(2); -c*y(2) + d*y(1)*y(2)];
[t, y] = ode45(rhs, [0, 100], y0);
plot(t, y, 'LineWidth', 3)
legend('Rabbits', 'Foxes')
Observing the periodicity of the solution, one may want to summarize it with a phase plot in addition to time-series plot:
figure()
plot(y(:,1), y(:,2), 'LineWidth', 3)
xlabel('Rabbits')
ylabel('Foxes')
