% Predator-Prey
% Consider the system of equations:
% \begin{array}{l}
% \frac{\textrm{dx}}{\textrm{dt}}=\textrm{kx}-\alpha \textrm{xy}\\
% \frac{\textrm{dy}}{\textrm{dt}}=-\textrm{cy}+\beta \textrm{xy}
% \end{array}
% with the inital conditions x(0)=1000 and y(0)=25. The built-in function ode45 will solve this system as follows

tr=[0 500];
initv=[1000 25];
[t,z]=ode45(@rkfunc, tr, initv);
plot(t,z(:,1), 'LineWidth', 2)
set(gca,'FontSize',20)
ylabel('x')
xlabel('time')
figure
plot(t,z(:,2), 'LineWidth', 2)
set(gca,'FontSize',20)
ylabel('y')
xlabel('time')
figure
plot(z(:,1),z(:,2), 'LineWidth', 2)
set(gca,'FontSize',20)
ylabel('x')
xlabel('y')

%
% The trange vector indicates that we will solve the equation on the range 0<t<500. The init vector sets the initial values of both x and y.
%
% As called, ode45 expects to find the function definition in a user-defined function which I've called rkfalling. For this problem, rkfalling is as follows:

function zout=rkfunc(t,z)
x=z(1);
y=z(2);
k=0.05;
c=0.018;
alpha=0.002;
beta=2e-5;
zout=zeros(2,1);
zout(1)=k*x-alpha*x*y;
zout(2)= -c*y+beta*x*y;
end

% Here z is a vector holding a value for both x and y. I've chosen x to be z(1) and y to be z(2). If this is the order of the variables in w, then rkfalling must return a vector which I've called r which holds values for dy/dt and dz/dt (in that order).
% © James Blanchard, University of Wisconsin-Madison, blanchard@engr.wisc.edu