# Phase Portraits using MATLAB
## Plot Solutions of a Planar ODE
In this section, we use ```pplane8.m``` to draw the direction field and some solutions of a two-dimensional system on the phase plane, even if the analytical solution cannot be found.

Some discussion of implementation errors may be found [here](https://www.mathworks.com/matlabcentral/answers/303035-question-re-pplane8-m).
1. Download the file [pplane8.m](https://www.mathworks.com/matlabcentral/fileexchange/61636-pplane).
2. Open pplane8.m and run it. (Or just type ```pplane8``` in the command window *when you are in the right folder.*) You should see the following window.

**Note.** If you want to setup the fontsize, copy the following code into the original file of ```pplane8.m```:

set(groot, 'DefaultUicontrolFontSize',18, 'DefaultAxesFontSize',18, 'DefaultTextFontSize',18)

In the following, we use the default example $$\begin{cases}
x'=2x-y+3(x^2-y^2)+2xy\\
y'=x-3y-3(x^2-y^2)+3xy.
\end{cases}$$
### Plot Some Solutions
1. Input the equations, choose the display window and the style of the direction field if you like. Then "Proceed" and you will see a new window with the direction field.
2. Click any point as you like to be an initial condition. It will plot the corresponding solution automatically. Of course you can draw more than one solutions. You can erase all solutions by clicking "Edit" and "Erase all solutions."
3. You can also draw a solution by inputing the initial condition using keyboard. First click "Solutions" and then "Keyboard input." You will see the following window and just input the initial condition as you like.
4. You may also plot the solution only in positive time. Click "Options" and "Solution direction" to set it.

### Equilibria and Closed Orbits
1. Click "Solutions" and then "Show nullclines." (This "Show" will be changed to "Hide.") The equilibria are on the intersection of two nullclines.
($x$-nullcline is the set $2x-y+3(x^2-y^2)+2xy=0$. Similarly for the $y$-nullcline.)
2. Click "Solutions" and then "Find an equilibrium point."
(An equilibrium is a constant solution.)
3. Click a point where seems to be near an equilibrium. Then it will plot a nearest equilibrium.
4. It simultaneously pops up a window (left) showing its Jacobian matrix at that equilibrium. You can display the linearization (right).
 

5. Similarly, click "Solutions" and "Find the nearly closed orbit." Then you can find a nearest closed trajectory (if it exists) when clicking a point.
### Varying Parameters
When we investigate a planar system with parameters, it is convenient to use the parameter row.

For example, consider the predator-prey system
$$\begin{cases}
x'=\gamma x(1-\dfrac{x}{K})-\alpha xy\\
y'=\beta xy-dy.
\end{cases}$$
with the parameters $\gamma,K,\alpha,\beta,d$. Then we can proceed as follows.
<p style="text-align: center;">
<img src="https://hackmd.io/_uploads/B1FMfI0gA.png" alt="pplanesetup" style="width:60%;">
</p>

Just type in other values of parameters to observe the discrepancy.
### Exercises
1. Consider the system
$$\begin{cases}
x'=-(x-y)(1-x-y)\\
y'=x(2+y).
\end{cases}$$Find all nullclines and equilibria, and plot them on the phase plane. Determine which equilibrium attracts other solutions and which does not.
2. Consider the system $$\begin{cases}
x'=\mu x+y-x(x^2+y^2)\\
y'=-x+\mu y-y(x^2+y^2).
\end{cases}$$
It is known that $(0,0)$ is the only equilibrium. Show numerically that all solutions tend to $(0,0)$ as $t\to\infty$ whenever $\mu<0$, and that all solutions tend to a closed orbit $x^2+y^2=\mu^2$ which is also known as a solution.
## Higher Dimensions
Find my codes in [this link](https://drive.google.com/drive/folders/1WizgC8OVL1YH0azcqI4O6ZTbypoH-xPw?usp=sharing).
### Systems of ODEs
For ODEs of higher dimensions
$${\bf x}'={\bf f}({\bf x}),~{\bf x}=(x_1,...,x_n),$$ we often only care about $x_i$-$t$ plot. For example, consider the SIR epidemic model
$$\begin{cases}
S'=\Lambda-\beta SI-dS\\
I'=\beta SI-\gamma I-dI\\
R'=\gamma I-dR.
\end{cases}$$
We first input the parameters and initial conditions. Then put them in vectors.

lambda = 1; d = 2e-3; beta = 4e-5; gamma = 2e-3;%parameters
pars = [lambda ; d; beta ; gamma];
S0 = 1000; I0 = 1; R0 = 0;%initial conditions
y0 = [S0 I0 R0];

Specify the time interval:

tspan = [0, 3000]; 

Now use an ODE solver (here we choose ```ode45```) in MATLAB to find the numerical solution.

[t,y] = ode45(@sir_rhs, tspan, y0, [], pars);

Here, ```sir_rhs``` is a function which carries the right-hand side of the ODE, which should be put in the last part of the code.

Next, we plot the solution. (We do not explain the technique of plotting here.)

figure(1);
plot(t,y(:,1),'b',t,y(:,2),'r',t,y(:,3),'k');
xlim(tspan);
ylabel('Population','FontSize',15);
xlabel('t','FontSize',15);
legend('$$S$$','$$I$$','$$R$$','Interpreter','Latex','Location','NorthEast','FontSize',15);

Go back to the function ```sir_rhs```. It has three inputs, namely, ```t,y,pars```, which are the time variable, components of the solution, and parameter vector.

function f = sir_rhs(t,y,pars)
f = zeros(3,1);
f(1) = pars(1)-pars(3)*y(2)*y(1)-pars(2)*y(1);
f(2) = pars(3)*y(2)*y(1)-pars(4)*y(2)-pars(2)*y(2);
f(3) = pars(4)*y(2)-pars(2)*y(3);
end 

Here we present the result.
<p style="text-align: center;">
<img src="https://hackmd.io/_uploads/ByiIYkWZA.jpg" alt="sir" style="width:60%;">
</p>

Of course, if the system is three-dimensional, we can plot the solution in the phase space by using ```plot3```.

figure(2);
plot3(y(:,1),y(:,2),y(:,3));
grid on;%if you don't want grid, just comment it.
xlabel('S','FontSize',15);
ylabel('I','FontSize',15);
zlabel('R','FontSize',15);

The result is presented below.
<p style="text-align: center;">
<img src="https://hackmd.io/_uploads/r1vNEdIZC.jpg" alt="vandepo" style="width:60%;">
</p>

We finally remark that ```ode45``` is usually used to to be our ODE solver. If the problem is *stiff* (namely, the solution is stable only in short time), it may not be an appropriate choice (see this [link](https://www.mathworks.com/help/matlab/math/choose-an-ode-solver.html)).

Appendix: full code (copy it directly)

lambda = 1; d = 2e-3; beta = 4e-5; gamma = 2e-3;%parameters
pars = [lambda ; d; beta ; gamma];
S0 = 1000; I0 = 1; R0 = 0;%initial conditions
y0 = [S0 I0 R0];
tspan = [0, 3000]; 
[t,y] = ode45(@sir_rhs, tspan, y0, [], pars);
figure(1);
plot(t,y(:,1),'b',t,y(:,2),'r',t,y(:,3),'k');
xlim(tspan);
ylabel('Population','FontSize',15);
xlabel('t','FontSize',15);
legend('$$S$$','$$I$$','$$R$$','Interpreter','Latex','Location','NorthEast','FontSize',15);
figure(2);
plot3(y(:,1),y(:,2),y(:,3));
grid on;%if you don't want grid, just comment it.
xlabel('S','FontSize',15);
ylabel('I','FontSize',15);
zlabel('R','FontSize',15);
function f = sir_rhs(t,y,pars)
f = zeros(3,1);
f(1) = pars(1)-pars(3)*y(2)*y(1)-pars(2)*y(1);
f(2) = pars(3)*y(2)*y(1)-pars(4)*y(2)-pars(2)*y(2);
f(3) = pars(4)*y(2)-pars(2)*y(3);
end

### Exercises
Consider the *Lorenz system*
$$\begin{cases}
x'=\sigma(y-z)\\
y'=\rho x-y-xz\\
z'=xy-\beta z.
\end{cases}$$
Fix $\sigma=10$ and $\beta=\dfrac{8}{3}.$ Do the following questions and explain what phenomena occur. Note that you need to use ```legend``` and color the curves to specify the components.
1. Fix $\rho=28$ and $(x(0),y(0),z(0))=(5,5,5)$. For $t\in[0,30]$, make the $x$-$t$, $y$-$t$ and $z$-$t$ plots in the same figure.
2. Fix $\rho=28$. For $t\in[0,30]$, make the $x$-$t$ plots under two different initial conditions $(x(0),y(0),z(0))=(5,5,5)$ and $(x(0),y(0),z(0))=(5.01,5,5)$ in the same figure.
3. Fix $\rho=21$. For $t\in[0,30]$, make the $x$-$t$, $y$-$t$ and $z$-$t$ plots under three different initial conditions $(x(0),y(0),z(0))=(3,8,0),(x(0),y(0),z(0))=(5,5,5)$ and $(x(0),y(0),z(0))=(5,5,10).$ (One figure for $(t,x(t)),(t,y(t)),(t,z(t))$ under one initial condition.)
4. Determine whether the solutions starting from a point near the origin change their behavior when moving $\rho$ in the interval $(\frac{1}{2},\frac{3}{2}).$ You should observe that a bifurcation occurs in this interval.
5. Plot some solutions (in different figures) in the phase space.

## Numerical Solutions of Higher-Order ODEs
Convert the ODE into a system of first-order ODEs (see [here](https://www.mathworks.com/help/matlab/ref/ode45.html#bu3uj8b)).
**Example.** The *van der Pol equation* is a second-order ODE $$y''-\mu(1-y^2)y'+y=0,$$where $\mu>0$. Let $x_1=y,x_2=y'$. Then
$$\begin{cases}
x_1'=x_2\\
x_2'=\mu(1-x_1^2)x_2-x_1.
\end{cases}$$

[t,y] = ode45(@vdp1,[0 20],[2; 0]);
plot(t,y(:,1),'-r',t,y(:,2),'-b')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution');
legend('y','dy/dt');
function dydt = vdp1(t,y)
%VDP1  Evaluate the van der Pol ODEs for mu = 1
dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];
end

We have
<p style="text-align: center;">
<img src="https://hackmd.io/_uploads/ByIVHwD-C.jpg" alt="vandepo" style="width:60%;">
</p>

Of course, you can ask ChatGPT to write down the codes for plotting the solutions, specifying your needs, including graph style, labels, etc.