How can I plot my solutions in my phase portait?

This model involves modeling Romeo and Juilet's love, however I am having trouble plotting my solutions. They are not showing in my direction field. I would like some help or insite on what's going on. I do notice there is an error that is affecting my plots.

% We start calculating the derivatives |y1'| and |y2'| for each point of the phase plane.
% We create a grid of points where we want to draw out plots:
y1 = linspace(-10,10,20);
y2 = linspace(-10,10,20);

%%
% |meshgrid| creates two matrices: one for all the uu-values of the grid, and
% one for all the vv-values of the grid. Then, we consider |x1| and |x2|
% matrices: |x1| will contain the value of |y1'| at each uu and vv position,
% while |x2| will contain the value of |y2'| at each uu and vv position of
% our grid.
[uu,vv] = meshgrid(y2,y1);
x1 = zeros(size(uu));
x2 = zeros(size(vv));

 a = 0.5; 
 b = 2; 
 c = 15; 
 d = 2; 
 e = 0.5; 
 f = 15; 
init_time=0;  
loveRJ = @(t,y) [a*y(2)*y(1) + b*y(1)*(1-(y(1)/c)); d*y(2)*y(1) + e*y(2)*(1-(y(2)/f))];

for i = 1:numel(uu)
    Yder = loveRJ(init_time,[uu(i); vv(i)]);
    x1(i) = Yder(1);
    x2(i) = Yder(2);
end

%%  
% Finally we compute a couple of solutions and plot them, along with the phase
% portrait, on the phase plane.
figure
quiver(gca,uu,vv,x1,x2,'r'); 
xlabel('Juliet Emotions');
ylabel('Romeo Emotions');
axis tight equal;

% Calculate and plot 1st Solution
tstart = 0;
tfinal = 50;
tspan = [tstart tfinal];

y0_1 = [2;-1]; % initial conditions
[t,y1] = ode23(@(t,y) loveRJ(t,y,a,b,c,d,e,f),tspan,y0_1);

Error using solution
Too many input arguments.

Error in solution (line 46)
[t,y1] = ode23(@(t,y) loveRJ(t,y,a,b,c,d,e,f),tspan,y0_1);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode23 (line 106)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

figure(gcf)
hold on
plot(y1(:,2),y1(:,1),'b')
plot(y1(1,2),y1(1,1),'mo',... % starting point
'MarkerEdgeColor','k',...
    'MarkerFaceColor',[.49 1 .63],...
    'MarkerSize',10)
plot(y1(end,2),y1(end,1),'bs',... % ending point
'MarkerEdgeColor','k',...
    'MarkerFaceColor',[.49 .63 1],...
    'MarkerSize',10)

% Calculate 2nd Solution
y0_2 = [4;1]; % initial conditions
[t,y2] = ode23(@(t,y)loveRJ(t,y,a,b,c,d,e,f),tspan,y0_2);
plot(y2(:,2),y2(:,1),'c')
plot(y2(1,2),y2(1,1),'ko',... % starting point
'MarkerEdgeColor','k',...
    'MarkerFaceColor',[.49 1 .63],...
    'MarkerSize',10)
plot(y2(end,2),y2(end,1),'bs',... % ending point
'MarkerEdgeColor','k',...
    'MarkerFaceColor',[.49 .63 1],...
    'MarkerSize',10)
hold off
% Calculate 3rd Solution
y0_3 = [5;2]; % initial conditions
[t,y3] = ode23(@(t,y) AgedloveRJ(t,y,a,b,c,d,e,f),tspan,y0_3);figure(gcf)
plot(y3(:,2),y3(:,1),'c')
plot(y3(1,2),y3(1,1),'ko',... % starting point
'MarkerEdgeColor','k',...
    'MarkerFaceColor',[.49 1 .63],...
    'MarkerSize',10)
plot(y3(end,2),y3(end,1),'bs',... % ending point
'MarkerEdgeColor','k',...
    'MarkerFaceColor',[.49 .63 1],...
    'MarkerSize',10)
hold off
title("Two solutions plotted on vector field")
title("Two solutions plotted on vector field")

Answer 1

[Steven Lord]

You may be interested in specifying the OutputFcn option in your ODE solver call. The odephas2 function included in MATLAB may do what you're looking for. See the orbitode example file for a demonstration of how to specify the OutputFcn as @odephas2.

[Jessica Arroyo]

what do you mean? Can you show me how to use it?