%% Direction Fields
%% Kurt Bryan and SIMIODE
% A script to illustrate how to draw a direction field for an ODE, and superimpose a solution curve.
%
% We'll use u'(t) = t*cos(u(t)) - sin(t) as an example. Define the right side t*cos(u)+sin(t) as a function of t and u:

f = @(t,u) t.*cos(u)-sin(t)

% Now construct the direciton field on the range 0 <= t <= 5, 0 <= u <= 5. First

[t,u] = meshgrid(0:0.25:5, 0:0.25:5);

% sets up a grid in the (t,u) plane on the given region with 0.25 spacing between grid points in each direciton. Then compute the slope "s" at each grid point with

s = f(t,u);

% Finally,have Matlab plot arrows at the appropriate grid points, with the appropriate slopes. We can scale the arrows to have an appealing length using a final argument "2" (you can choose it as you wish):

plot1 = quiver(t,u,ones(size(s)), s, 2);
axis tight

% To superimpose solution curves, solve the ODE numerically. For example, to show a solution curve with initial data u(1) = 2, call Matlab's solver

inittime = 1.0; ic = 2.0; finalT = 5.0;
[tm,um] = ode45(f,[inittime,finalT], ic); %Matlab's solver
hold on;
plot2 = plot(tm,um,'-r');
hold off

% Alternatively, use the supplied subroutine "draw_dirfield" in the form draw_dirfield(f, range ,ics) where "range" is of the form [lowt, hight, lowu, highu] and ics is an optional n x 2 matrix of initial conditions with kth row tk, uk (where u(tk) = uk.)

ics = [0 1;1 2;1 4];
draw_dirfield(f,[0,5,0,5],ics)