How to set a point which coordinates depends on two functions

Basically I'm resolving a system of two differential equations with a R-K of order 4 and the point of equilibrium y0 from which the solutions begins isn't fixed. In fact y0 depends from two functions appearing in the system. So the point y0 stays on a curve defined from these 2 functions which I've mentioned before. What I'm asking for is if there's a way of plotting the solutions with y0 running on this curve.

My codes are:

function [t,y]=RK(fun,tspan,y0,n,A,b,c)
h=(tspan(2)-tspan(1))/n;
y(:,1)=y0;
t(1)=tspan(1);
s=length(b); 
for j=1:n 
    k(:,1)=feval(fun,t(j),y(:,j));
    for i=2:s
        k(:,i)=feval(fun,t(j)+c(i)*h,y(:,j)+h*k(:,1:i-1)*A(i,1:i-1)');
    end
    y(:,j+1)=y(:,j)+h*k*b;
    t(j+1)=t(j)+h;
end
end

the first one

function dy=pred_prey(t,y)
a=1/3;
d=2/3;
r=(2/3)*(1/(exp(-0.5*t)+1)-0.5)*(1+(sin(t)/t));
mu=13/20-(3/5)*exp(-(1/(2*t)));
dy(1)=y(1)*r-a*y(1)*y(2);
dy(2)=-mu*y(2)+d*y(1)*y(2);
dy=dy';
end

second one

function []=driver_pred_prey()

close all;
tspan=[1/100 50];
n=5000;
fun='pred_prey';
t=[0.01:0.01:50];
r=(2/3)*(1./(exp(-0.5*t)+1)-0.5).*(1+(sin(t)./t));
mu=13/20-(3/5)*exp(-(3./t));
y0=[0.83705 0.45454]';

A=[0 0 0 0; 1/2 0 0 0 ; 0 1/2 0 0; 0 0 1 0];
b=[1/6 1/3 1/3 1/6]';
c=[0 1/2 1/2 1]';

[t,y]=RK(fun, tspan, y0, n, A, b, c);

figure (1)
plot(t, y(1,:), 'k-'), hold
plot(t, y(2,:), 'k-.') 
grid on
legend('resources','population')

figure (2)
plot(y(1,:),y(2,:))
title('orbits')

end

It runs obviously. I've got the y0 that you find here from plotting separatley the curve it stays on, and then getting its coordinates.