First download the file `vectfield.m`

.

Here we consider the following example of an autonomous system:

y_{1}' = 4y_{1}+y_{2}+y_{1}y_{2}

y_{2}' =y_{1}+ 4y_{2}+y_{2}^{2}

**First define the right hand side function
**

of the differential equation as**f**

`f = @(t,y) [4*y(1)+y(2)+y(1)*y(2) ; y(1)+4*y(2)+y(2)^2]`

Note that you **must **use

, despite the fact that the
functions do not depend on **@(t,y)** ...*t*.

**To plot the vector field** for *y*_{1}
between -6 and 5 and *y*_{2} between -6 and 1 type

We now want to plot the trajectory through a point (a1,a2), say (a1,a2)=(1,-3):

`vectfield(f,-6:.5:5,-6:.5:1)`

a1=1; a2=-3

**To plot a trajectory in the phase
plane** starting at a point (`a1`

, `a2`

) at
time *t*=0 for *increasing* values of *t* going from 0 to
4 type

[ts,ys] = ode45(f,[0,4],[a1;a2]); plot(ys(:,1),ys(:,2))

To plot the trajectory in the phase plane starting at the point
(`a1`

, `a2`

) at time *t*=0 for
*decreasing* values of *t* going from 0 to -4 type

[ts,ys] = ode45(f,[0,-4],[a1;a2]); plot(ys(:,1),ys(:,2))

To get an idea of the behavior of the ODE, we **plot the direction
field and several trajectories together**: We can choose e.g. the 63
starting points (`a1`

, `a2`

) with
`a1`

=-7,-5,-3,-1,1,3,5 and `a2`

=-7,-6,-5,-4,-3,-2,-1,0,1. From each
of those points we plot the trajectory for *t* going from 0 to 4 and
for *t* going from 0 to -4:

vectfield(f,-6:.5:5,-6:.5:1) hold on for a1=-7:2:5 for a2=-7:1 [ts,ys] = ode45(f,[0,4],[a1;a2]); plot(ys(:,1),ys(:,2)) [ts,ys] = ode45(f,[0,-4],[a1;a2]); plot(ys(:,1),ys(:,2)) end end hold off

These commands will generate a number of warning messages (since some
solutions don't exist for *t* in the whole interval [-4,4]) which you
can ignore. You can use the command

warning('off','MATLAB:ode45:IntegrationTolNotMet')at the beginning of your m-file to suppress these warnings. Then your code should run without any warnings or errors. From the picture it appears that there are three critical points:

- a stable node at (-5,-5)
- an unstable node at (0,0)
- a saddle point at (3,-3)

Recall that we defined the right hand side function

of the differential equation as**f**

`f = @(t,y) [4*y(1)+y(2)+y(1)*y(2) ; y(1)+4*y(2)+y(2)^2]`

We define symbolic variables

syms t y1 y2

and define `F`

by using `f`

with the vector
`[y1;y2]`

F = f(t,[y1;y2])

which gives the symbolic vector

[ 4*y1+y2+y1*y2] [ y1+4*y2+y2^2]

**To find critical points** where both components
`F(1)`

and `F(2)`

of `F`

are zero we use

:**solve**

[y1s,y2s] = solve(F(1),F(2),y1,y2)

This returns

y1s = [ 0] [ -5] [ 3] y2s = [ 0] [ -5] [ -3]

Therefore we have found three critical points: (0,0), (-5,-5) and (3,-3). Your version of Matlab may give the answer in a different order.

**To find the Jacobian matrix** we type

`jac = jacobian(F,[y1;y2])`

This returns the matrix

[ 4+y2, 1+y1] [ 1, 4+2*y2]

As an example we will study the second critical point with coordinates
`y1s(2)`

=-5, `y2s(2)`

=-5.

**To evaluate the Jacobian at the critical point** we
substitute `y1s(2)`

and `y2s(2)`

for `y1`

and
`y2`

:

`A = subs(jac,{y1,y2},{y1s(2),y2s(2)})`

which gives the matrix

[ -1, -4] [ 1, -6]

**To find the eigenvectors and eigenvalues of this matrix**
type

`[eigvect,eigval] = eig(A)`

which gives

eigvect = [ 1, 4] [ 1, 1] eigval = [ -5, 0] [ 0, -2]

This means that the eigenvalues are -5 and -2. Using the table for critical points we see that this critical
point is a *stable node*. The eigenvector for the first eigenvalue -5
is given by the first column of the matrix `eigvect`

, it has
components 1, 1. The eigenvector for the second eigenvalue -5 is given by the
second column of the matrix `eigvect`

, it has components 4, 1.
Using the eigenvectors we can sketch by hand the behavior of the trajectories
close to the critical point:

Compare this picture with the trajectories near (-5,-5) in the picture above.

You should try to analyze the other two critical two points by yourself in the same way.