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 or qualitative information about general behavior of the solution curves (call trajectories or streamlines) from the differential equation itself without actual solution formula.
A direction field or a slope field for a first order differential equation \( {\text d}y / {\text d}x = f(x,y) , \) is a field of short either straight line segments or arrows of slope f(x,y) drawn through each point (x,y) in some chosen grid of points in the (x,y) plane.

Direction fields could be visualized by plotting a list of vectors that are tangent to solution curves. It is as if we had a traffic policeman stationed at each point (x,y) on the plane and directed the traffic to flow along the direction specified by the slope function f(x,y). From this point of view, a solution is a curve that obeys the law at each of its points. Correspondingly, 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 𝑎 ≤ xb, cyd, type the following sequence of commands:


[x, y] = meshgrid(-2:0.2:3, -1:0.2:2); 
s = 1 - x.*y.^2;    % for slope function f = 1 - x y^2
quiver(x,y, ones(size(s)),s), axis tight
The first command sets 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 likr this:

picture??????????????

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:


[x, y] = meshgrid(-2:0.2:3, -1:0.2:2); 
s = 1 - x.*y.^2;    % for slope function f = 1 - x y^2
L = sqrt(1 + s.^2);)
quiver(x,y, 1./L,s./L, 0.5), axis tight
xlabel 'x', ylabel 'y' 
title 'Direction field for dy/dx = 1-xy^2'
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

picture??????????????

Here is another matlab code using trigonometric functions. We split the slope into x and y components and normalize it using arctangent, sine, and cosine, allowing us to plot vectors of magnitude 1 for our direction field. Finally, we plot the separatrix, which needs to be derived by hand and inputted into matlab.

% Meshgrid sets up the x-coordinate range and y-coordinate range
[x,y] = meshgrid(0:0.2:3, -1:0.2:3);

% Next, set up matrix of slope values from the differential equation
% Be sure to yse .* instead of * for element-wise matrix multiplication
% (as opposed to matrix-matrix multiplication)
diff_val = 5.*y - 3.*x;

% Calculate the angle theta between the slope and the horizontal using the
% arctangent function
theta = atan(diff_val);

% Using theta, you can split the slope into x and y components (note that
% the vector [u,v] is now a unit vector
u = cos(theta);
v = sin(theta);

% Initialize a figure
figure

% Use the quiver function to plot slope vectors [u,v] for the ranges [x,y]
% as defined in the meshgrid function
quiver(x,y,u,v)

% To also graph the separatrix, use the hold on command, which will add
% additional plots to the same figure. Note that you must find the
% separatrix by hand and then enter its equation into MATLAB
hold on
fplot(@(a) (1/25)*(3+15*a),[0 3],'g')

% Below is some code for labeling our graph and adding a title
xlabel('x'); ylabel('y'); title('Slope field for dy/dx = 5y - 3x');
hold off

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.


[x, y] = meshgrid(-3:.3:3, -2:.3:2);

% defines ODE
dy = x + sin(y); 
dx = ones(size(dy));

% normalizes vectors
norm_factor = sqrt(dx.^2 + dy.^2);
dy_norm = dy./norm_factor;
dx_norm = dx./norm_factor;

quiver(x, y, dx_norm, dy_norm);
Another code with numerical solution

%  portrait.m 
[x,y] = meshgrid(-3:0.2:3,-3:0.2:3);
n = size(x);
u = ones(n);
v = 1 - x.*y.^2; 

for 1 = 1:numel(x)
rac = 2*sqrt(u(i)^2 + v(i)^2); 
u(i) = u(i)/rac; 
v(i) = v(i)/rac;
end

xspan = [-2,3];
x0 = -2; 
[xx,yy] = ode45(@(x,y) 1-x*y^2, xspan, x0);
figure
hold on 
quiver(x,y,u,v)
plot(xx,yy,'-')

The best way to plot direction fileds 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 fileds 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 . \)

% plot direction field for first order ODE y' = f(t,y)
function dirfield
              
f = @(t,y) -y^2+t % f = @(independent_var, dependent_var)
    % mathematical expression y in terms of t

figure;             
dirplotter(f, -1:.2:3, -2:.2:2) % dirfield(f, t_min:t_spacing:t_max, y_min:y_spacing:y_min)
    % using t-values from t_min to t_max with spacing of t_spacing
    % using y-values from y_min to t_max with spacing of y_spacing
end

function dirplotter(f,tval,yval) %  plots direction field using aboved defined (f,tval,yval)
[tm,ym]=meshgrid(tval,yval); % creates a grid of points stored in matrices in the (t,y)-plane
  dt = tval(2) - tval(1); % determines t_spacing and y_spacing based on tval
  dy = yval(2) - yval(1); % determines t_spacing and y_spacing based on yval
  fv = vectorize(f); % prevents matrix algebra; if f is a function handle fv is a string
  if isa(f,'function_handle') % if f is function handle
  fv = eval(fv); % turn fv into a function handle from a string 
  end
  yp=feval(fv,tm,ym); % evaluates differential equation fv from tm and ym from meshgrid
  s = 1./max(1/dt,abs(yp)./dy)*0.35;
  h = ishold; % define h as ishold which returns 1 if hold is on, and 0 if it is off
  quiver(tval,yval,s,s.*yp,0,'.r'); hold on; % plots vectors at each point w direction based on righthand side of differential equation
  quiver(tval,yval,-s,-s.*yp,0,'.r');
  if h % h above defined as ishold
  hold on % retain current plot 
  else % refers to all not h 
  hold off % resets axes properties to default
  end
  axis([tval(1)-dt/2,tval(end)+dt/2,yval(1)-dy/2,yval(end)+dy/2]) % sets the scaling for axes on current graph to ends of tval and yval   
end    

% 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. Babolian, E. and Javadi, Sh., “New method for calculating Adomian polynomials”, Applied Mathematics and Computation, 2004, Volume 153, Issue 1, 25 May 2004, pages 253--259. https://doi.org/10.1016/S0096-3003(03)00629-5