# How to generate phase plots for the solutions of this ODE system?

tmax = 500;
\[Beta]a = 0.4;
\[Beta]i = 0.3;
b = 0.55;
\[Nu] = 0.3;
\[Mu] = 0.01;
\[Rho] = 0.1;
\[Xi] = 0.001;
\[Delta]a = 0.15;
\[Delta]i = 0.41;
\[Sigma] = 0.8;
\[Tau] = 0.6;
\[Alpha] = 0.01;
SAExR = NDSolveValue[{
S'[t] ==
b*(1 - \[Nu]) - (\[Beta]a *A[t] + \[Beta]i *Ex[t])*
S[t] - (\[Mu] + \[Rho])*S[t] + \[Xi] *R[t],
A'[t] == (\[Beta]a* A[t] + \[Beta]i* Ex[t])*
S[t] - (\[Mu] + \[Delta]a + \[Sigma])*A[t],
Ex'[t] == \[Sigma]*A[t] - (\[Mu] + \[Delta]i + \[Tau] + \[Alpha])*
Ex[t],
R'[t] ==
b*\[Nu] + \[Rho]*S[t] + \[Delta]a*A[t] + (\[Delta]i + \[Tau])*
Ex[t] - (\[Xi] + \[Mu])*R[t],
S[0] == 30,
A[0] == 10,
Ex[0] == 10,
R[0] == 0},
{S, A, Ex, R},
{t, 0, tmax}];
{f1, f2, f3, f4} = SAExR;

st = Style[#, 15, Black] &;

Plot[{f1[t], f2[t], f3[t], f4[t]}, {t, 0, 50},
PlotStyle -> {Blue, Orange, Red, Green}, Frame -> True,
FrameLabel -> {Style["Time", 20, Black],
Style["Density", 20, Black]},
PlotLegends ->
Placed[LineLegend[{Blue, Orange, Red, Green}, {"S(t)", "A(t)",
"I(t)", "R(t)"}, LegendFunction -> Framed], {0.85, 0.3}],
ImageSize -> 500, FrameTicksStyle -> 18, FrameStyle -> Black]
Plot[{f1[t], f2[t], f3[t]}, {t, 0, 50},
PlotStyle -> {Blue, Orange, Red}, Frame -> True,
FrameLabel -> {Style["Time", 20, Black],
Style["Density", 20, Black]},
PlotLegends ->
Placed[LineLegend[{Blue, Orange, Red}, {"S(t)", "A(t)", "I(t)"},
LegendFunction -> Framed], {0.9, 0.8}], ImageSize -> 500,
FrameTicksStyle -> 18, FrameStyle -> Black]
Plot[{f4[t]}, {t, 0, 500}, PlotStyle -> {Green}, Frame -> True,
FrameLabel -> {Style["Time", 20, Black],
Style["Density", 20, Black]},
PlotLegends ->
Placed[LineLegend[{ Green}, {"R(t)"},
LegendFunction -> Framed], {0.9, 0.75}], ImageSize -> 500,
FrameTicksStyle -> 18, FrameStyle -> Black]


How would I generate the phase plots for this system? For the example (R,S) phase plot?

Edit: Can someone create something like this but with the same plotting structure(same axis numbering, in black, same font size etc) as in the numerical simulations for the ODE system.

Normally, a phase plot displays a function plotted against its derivative, {f,f'}. You ask for a "phase plot" of R,S, but are they derivatives of each other? Or are you asking for a parametric plot of R[t] vs. S[t]? What about Ex and A? Or {f,f'}? I'll provide you with a set of each.

I used NDSolve in place of NDSolveValue with the same constants:

SAExR = NDSolve[{
S'[t] == b*(1 - \[Nu]) - (\[Beta]a *A[t] + \[Beta]i *Ex[t])* S[t] - (\[Mu] + \[Rho])*S[t] + \[Xi] *R[t],
A'[t] == (\[Beta]a* A[t] + \[Beta]i* Ex[t])*S[t] - (\[Mu] + \[Delta]a + \[Sigma])*A[t],
Ex'[t] == \[Sigma]*A[t] - (\[Mu] + \[Delta]i + \[Tau] + \[Alpha])*Ex[t],
R'[t] == b*\[Nu] + \[Rho]*S[t] + \[Delta]a*A[t] + (\[Delta]i + \[Tau])*Ex[t] - (\[Xi] + \[Mu])*R[t],
S[0] == 30, A[0] == 10, Ex[0] == 10, R[0] == 0},
{S, A, Ex, R},{t, 0, tmax}];


Some of the options in your original plot are unnecessary for a ParametricPlot, so here's what's left after removing stuff like PlotLegends, etc.

GraphicsRow[{ParametricPlot[{S[t],R[t]}/.SAExR, {t, 0, 50},
PlotStyle ->Blue, Frame -> True,
FrameLabel -> {Style["S[t]", 20, Black],
Style["R[t]", 20, Black]},
ImageSize -> 500, FrameTicksStyle -> 18, FrameStyle -> Black],
ParametricPlot[{Ex[t],A[t]}/.SAExR, {t, 0, 50},
PlotStyle -> Orange, Frame -> True,
FrameLabel -> {Style["Ex[t]", 20, Black],
Style["A[t]", 20, Black]},
ImageSize -> 500,
FrameTicksStyle -> 18, FrameStyle -> Black]}]


ParametricPlot can phase plots (derivative of function vs function) as well. I modified your image size and aspect ratio, since the values you chose for Plot made them hard to copy/paste from my laptop.

 ParametricPlot[{S[t],S'[t]}/.SAExR, {t, 0, 50},
PlotStyle ->Blue, Frame -> True,
FrameLabel -> {Style["S[t]", 20, Black],
Style["S'[t]", 20, Black]},
ImageSize -> 350, FrameTicksStyle -> 18, FrameStyle -> Black,AspectRatio->2]


ParametricPlot[{R[t],R'[t]}/.SAExR, {t, 0, 50},
PlotStyle -> Orange, Frame -> True,
FrameLabel -> {Style["R[t]", 20, Black],
Style["R'[t]", 20, Black]},
ImageSize -> 350, FrameTicksStyle -> 18, FrameStyle -> Black,AspectRatio->2]


 ParametricPlot[{Ex[t],Ex'[t]}/.SAExR, {t, 0, 50},
PlotStyle ->Red, Frame -> True,
FrameLabel -> {Style["Ex[t]", 20, Black],
Style["Ex'[t]", 20, Black]},
ImageSize -> 350, FrameTicksStyle -> 18, FrameStyle -> Black,AspectRatio->2]


ParametricPlot[{A[t],A'[t]}/.SAExR, {t, 0, 50},
PlotStyle -> Green, Frame -> True,
FrameLabel -> {Style["A[t]", 20, Black],
Style["A'[t]", 20, Black]},
ImageSize -> 350, FrameTicksStyle -> 18, FrameStyle -> Black,AspectRatio->2]