For any differential equation in normal form

\[ y' = f(x,y) \qquad \mbox{or} \qquad \frac{{\text d}y}{{\text d}x} = f(x,y) \]
where f(x,y) is a well defined in some domain slope function, it is possible to obtain a graphical information about general behavior of the solution curves (call trajectories) from the differential equation itself without actual solution formula.

In the previous section, we show that direction fields or slope fields are very important features of differential equations because they provide a qualitative behavior of solutions. However, more precise information results from including in the plot some typical solution curves or trajectories.

A plot that shows representative sample of trajectories for a given first order differential equation is called phase portrait.
This section shows how to include sample trajectories into tangent field to obtain a phase portrait for a given differential equation.

The main matlab command for plotting direction fields is quiver, used in conjuction with meshgrid. To plot the slope field of a differential equation \( y' = f(x,y) \) on the rectangle 𝑎 ≤ x ≤ b, c ≤ y ≤ d, type the following sequence of commands:

The first command sets up a 26 by 16 grid of uniformly spaced points in the rectangular domain [-2,3] × [-1,2]. The second command evaluates the right-hand side of the differential equation \( y' = 1- x\,y^2 \) at the grid points, thereby producing the slopes, s = f(x,y), at these points. The vectors (1,s) have the desired slopes but have different lengths at the different grid points. The quiver command, used for plotting vector fields, requires four inputs: the array x of x-values, the array y of y-values, and arrays consisting of the two components of the direction field vectors. Since all of these arrays must have the same size, the code ones(size(s)) conveniently creates an array of ones of the same size as s. In the last command, axis tight eliminates white space at the edges of the direction field. The result looks like this:

While this picture is correct, it is somewhat hard to read because of the fact that many of the vectors are quite small. You can get a better picture by rescaling the arrows so that they do not vary in magnitude---for example, by dividing each vector (1,s) by its length \( \| (1,s)\| = \sqrt{1 + s^2} . \) In order to achieve it in matlab, we modify the preceeding sequence of commands to:

      The fifth entry, 0.5, in the quiver command reduces the length of the vectors by half and prevents the arrow heads from swallowing up the tails of nearby vectors. The result is presented at the left.

The picture strongly suggests that some solutions approach a specific curve, but others tend to minus infinity. This specific curve is called the separatrix.


Unless the special purposed field program is used, plotting slope fields in matlab requiries a bit of work. The following program utilizes the built-in command quiver(x,y,dx,dy) that plots an n × n array of vectors with x- and y- components specified by the n×n matrices dx and dy, based at the xy-points in the plane whose x- and y-coordinates are specified by the n×n matrices x and y.

Another Octave code:
Another code with numerical solution

We use Chebfun to plot phase portrait:
??????
Error using chebop/solveivp (line 240) The number of initial/final conditions does not match the differential order(s) appearing in the problem. Error in \ (line 49) [varargout{1:nargout}] = solveivp(N, rhs, pref, varargin{:});
    

N = chebop(0,100);
N.op = @(t,u) diff(u, 1) - 3*(1-u^2) + t.*u;

% Then we call quiver with N as an argument
quiver(N, [-2.75 2.75 -5.5 5.5],'xpts', 40, 'ypts', 40, 'scale', .5, ...
        'normalize', true)
hold on
for init = 0.2:0.4:0.2
    N.lbc = init;
    u = N\0;
    plot(1, u)
end
title('Phase portrait of the first-order ODE','interpreter','latex', ...
    'fontsize',11)
xlabel('$u$','interpreter','latex', 'fontsize',10)
ylabel('$u''$','interpreter','latex', 'fontsize',10)
hold off

To plot direction fileds along with some solutions, use the command streamline 
 % define arrays x, y, u, and v:
[x,y] = meshgrid(0:0.1:1,0:0.1:1);
u = x;
v = -y;
% Create a quiver plot of the data. Plot streamlines that start at different points along the line y=1.
figure
quiver(x,y,u,v)

startx = 0.1:0.1:1;
starty = ones(size(startx));
streamline(x,y,u,v,startx,starty)

The best way to plot direction fields is to use existing m-files, credited to John Polking from Rice University (http://math.rice.edu/~polking/). He is the author of two special matlab routines: dfield8 plots direction fields for single, first order ordinary differential equations, and allows the user to plot solution curves; pplane8 plots vector fields for planar autonomous systems. It allows the user to plot solution curves in the phase plane, and it also enables a variety of time plots of the solution. It will also find equilibrium points and plot separatrices. The latest versions of dfield8 and pplane8 m-functions are not compatible with the latest matlab version. You can grab its modification that works from dfield.m

Since matlab has no friendly subroutine to plot direction fields for ODEs, we present several codes that allow to plot these fields directly. Many others can found on the Internet. The following Octave code shows how to plot the direction field for the linear differential equation \( y' = 5\,y - 3\,x . \)

The following code presents the function version to plot direction field for first order ODE y' = f(t,y)

Another function: 

% Simple direction field plotter
clear;
[T Y]=meshgrid(-2:0.2:2,-2:0.2:2);
dY=cos(3*T)+Y./T;
dT=ones(size(dY));
scale=sqrt(1+dY.^2);
figure;
quiver(T, Y, dT./scale, dY./scale,.5,'Color',[1 0 0]);
xlabel('t');
ylabel('y');
title('Direction Field for dy/dt');
axis tight;
hold on;

% change color
% do multiple examples, x.^2-y.^2; use T, etc

% Simple direction field plotter
clear;
[T Y]=meshgrid(-2:0.2:2,-2:0.2:2);
dY=T.^2-Y.^2;
dT=ones(size(dY));
scale=sqrt(1+dY.^2);
figure;
quiver(T, Y, dT./scale, dY./scale,.5,'Color',[0 1 0]);
xlabel('t');
ylabel('y');
title('Direction Field for dy/dt');
axis tight;
hold on;


			Complete
  1. Trefethen, L.N., Birkisson, A., Driscoll, T.A., Exploring ODEs, SIAM, Philadelphia, 2018.