function predprey

x0 = input('Enter the initial x-value: ');
y0 = input('Enter the initial y-value: ');
t0 = input('Enter the initial time: ');
tfinal = input('Enter the end time: ');
N = input('Enter the number of intervals: ');

h = (tfinal - t0)/N;
t = zeros(N+1,1);
x = zeros(N+1,1);
y = zeros(N+1,1);
t(1) = t0;
x(1) = x0;
y(1) = y0;
for i = 1:N
    k1x = evalf(t(i),x(i), y(i));
    k1y = evalg(t(i),x(i),y(i));
        
    k2x = evalf(t(i)+0.5*h,x(i)+0.5*h*k1x,y(i)+0.5*h*k1y);
    k2y = evalg(t(i)+0.5*h,x(i)+0.5*h*k1x,y(i)+0.5*h*k1y);
    
    k3x = evalf(t(i)+0.5*h,x(i)+0.5*h*k2x,y(i)+0.5*h*k2y);
    k3y = evalg(t(i)+0.5*h,x(i)+0.5*h*k2x,y(i)+0.5*h*k2y);
    
    k4x = evalf(t(i)+h,x(i)+h*k3x,y(i)+h*k3y);
    k4y = evalg(t(i)+h,x(i)+h*k3x,y(i)+h*k3y);
    
    x(i+1) = x(i) + h/6*(k1x + 2*k2x + 2*k3x + k4x);
    y(i+1) = y(i) + h/6*(k1y + 2*k2y + 2*k3y + k4y);    
    t(i+1) = t(i) + h;
end
%theoutput = [t x y]
plot(t,x,t,y,'--')
xlabel('time')
ylabel('Population')
title('Predator Prey Problem')
legend('Prey', 'Predator')
% end function predprey

function f = evalf(t,x,y)
f = 0.2*x - 0.005*x*y;
% end function evalf

function g = evalg(t,x,y)
g = - 0.5*y + 0.01*x*y;
% end function evalg