For any differential equation in normal form

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 𝑎 ≤ x ≤ b, c ≤ y ≤ d, 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

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'

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

% 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);

%  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

figure(1)
plot(x, y, 'Color', [1 0 0]) %blue line
hold on
plot(x, z, 'Color', [0 1 0]) %green line