%%%%This program solves Romio and Juliet model, plots the phase plane, and
%%%%plots the phase portrait
clear all
close all
global a b c d alpha beta
a = 0;  %%The rate of change of Romeo's love WRT romeo's love
b = 1;  %%The rate of change of Romeo's love WRT juliet love
c = -1; %%The rate of change of Julie't love WRT romeo's love
d = 0;  %%The rate of change of Juliet's love WRT julie't love

[T,Y] = ode45(@romio,[0 20],[1 1]);

subplot(1,2,1)
plot(Y(:,1),Y(:,2),'r','linewidth',2);
xlabel('Romeo"s love for Juliet')
ylabel('Juliet"s love for Romeio')
title('Linear Model of Love Affair')
hold on
%%%Plotting the romio and juliet phase portrait
%%%First, we calculate a mesh grid that composes the Romio-Juliet space
[R,J] = meshgrid(-2:0.25:2,-2:0.25:2);

%%%%Next, we calculate the "Velocity" of each point in Romio-Juliet space
%%%%of both Romio and Juliet
dRdt = a.*R + b.*J;
dJdt = c.*R + d.*R;

%%%%We use the matlab function 'quiver' to plot the total phase
%%%%portrait
subplot(1,2,1)
quiver(R,J,dRdt,dJdt)
hold off

%%%Non-linear love affair
%%%Now, we solve the new system described in the .m file 'romiononlinear,m'
%%%%The parameters for the nonlinear model are:
alpha = 1;
beta = 1;

%%%%First, we solve and plot the phase plane
[TNL,YNL] = ode45(@romiononlinear,[0 8],[.01 .5]);
subplot(1,2,2)
plot(YNL(:,1),YNL(:,2),'r','linewidth',2);
xlabel('Romeo"s love for Juliet')
ylabel('Juliet"s love for Romeio')
title('Non-Linear Model of Love Affair')
hold on

%%%%Now, we solve and plot the phase portrait
[R,J] = meshgrid(0:2:30,0:2:30);
dRnldt = alpha.*R;
dJnldt = -beta.*(R.^2 - 1).*J - R;
subplot(1,2,2)
quiver(R,J,dRnldt,dJnldt,3)
hold off
