# Using Matlab for Autonomous Systems

First download the file `vectfield.m` .

Here we consider the following example of an autonomous system:

y1' = 4y1 + y2 + y1 y2

y2' = y1 + 4 y2 + y22

## Plotting the vector field and trajectories

First define the right hand side function `f` of the differential equation as

```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 `@(t,y) ...`, despite the fact that the functions do not depend on t.

To plot the vector field for y1 between -6 and 5 and y2 between -6 and 1 type

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

We now want to plot the trajectory through a point (a1,a2), say (a1,a2)=(1,-3):
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)

## Investigation of critical points

Recall that we defined the right hand side function `f` of the differential equation as

```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.