state L=0, P=0, Z=0

param alpha1 = 3
param alpha2 = 1
param alpha3 = 0.1

param beta1 = 1
param beta2 = 5
param beta3 = 10

param gamma = 1
param delta = 1
param A_L = 2
param A_P = -1

L'= -alpha1*L + beta1*(P*(1 - (P/gamma)^2) + A_P)
P'= -alpha2*P + beta2*(L + A_L/(1 + delta*Z))
Z'= -alpha3*Z + beta3*P

//print t, L, P, Z
plot title "Rinaldi Love Dynamics Model" legend grid xlabel="Time" style="lines" t, L title="L", P title="P", Z title="Z"
plot title "L-P Phase Plane" grid xlabel="L" ylabel="P" L, P

option tstart=0, tstop=20, outstep=-0.1