function prob_2_29(fun,window,npts)
% Form a direction field and plot solutions for a scalar ODE.
% FUN is an expression in the form of a string for evaluating the 
% differential equation y' = f(t,y). The independent variable must 
% be t and the dependent variable must be y. 
% WINDOW = [tl tr yb yt] says that the plot is to take place in a 
% window of tl <= t <= tr and yb <= y <= yt. If WINDOW is not 
% supplied, it is given the default value of [-10 10 -10 10].
% NPTS = [nt ny] says that the direction field is to have nt equally
% spaced points in t and ny equally spaced points in y. If NPTS is
% a scalar, both nt and ny are given that value. If NPTS is not 
% supplied, it is given the default value of 20.
% It is assumed that the differential equation is moderately stable.
% The command
%>> dfs('cos(t)*y',[0 12 -6 6])
% forms a direction field for y' = cos(t)*y and readies the figure
% for input of initial conditions for a solution. Fig. 1.4 of the
% book could be reproduced in this way.
% Clicking at a point within the window causes the solution through
% the point to be computed and plotted up to the edge of the window.
% To cope with singular points, the computation is terminated where a
% small step size is needed to resolve a sharp change in the solution. 
global wB wT hmin f
% Process input data.
if nargin < 3
    nt = 20;
    ny = 20;
elseif length(npts) == 1
    nt = npts(1);
    ny = npts(1);
else
    nt = npts(1);
    ny = npts(2);
end
if nargin < 2
    window = [-10 10 -10 10];
end
wL = window(1);
wR = window(2);
wB = window(3);
wT = window(4);
if wR <= wL
    error('Must have window(1) < window(2).')
elseif wT <= wB
    error('Must have window(3) < window(4).')
end
hmin = 1e-4*(wR - wL);
f = vectorize(inline(fun,'t','y'));
% Set up grid of points for arrows of direction field:
tpts = linspace(wL,wR,nt);
ypts = linspace(wB,wT,ny);
[T,Y] = meshgrid(tpts,ypts);
% Evaluate the slopes and then normalize the vectors (1,S(i,j)). 
% Note that f has been vectorized so as to do this efficiently.
S = f(T,Y);
L = sqrt(1 + S.^2);
NS = S ./ L;
NL = 1 ./ L;
% The arrows are long, so shorten them by half.
quiver(T,Y,NL,NS,1/2,'b')
axis tight
title('Click at a point inside the window to get a solution.')
xlabel('Click at a point outside the window to quit.')
% Use an event function to terminate the integration on going
% out the bottom or top of the window. Use an output function
% to trap singular points by terminating on a minimum step size.
options = odeset('OutputFcn',@outfcn,...
                 'RelTol',1e-4,'AbsTol',1e-7*max(abs(wB),abs(wT)));
hold on
while 1
    % Get the initial condition from the plot.
    [t0,y0] = ginput(1);
    % Quit when the initial condition is outside the window.
    if (t0 <= wL) | (wR <= t0) | (y0 <= wB) | (wT <= y0)
        break;
    end
    [t,y] = ode45(@F,[t0,wR],y0,options);
    plot(t,y,'r')
    [t,y] = ode45(@F,[t0,wL],y0,options);
    plot(t,y,'r')
end
hold off
%===================================================
function yp = F(t,y)
global f
yp = f(t,y);

% $$$ function [value,isterminal,direction] = events(t,y)
% $$$ global wB wT
% $$$ value = [wB; wT] - y;
% $$$ isterminal = [1; 1];
% $$$ direction = [0; 0];


function status = outfcn(t,y,flag)
global hmin wB wT 
persistent previoust
status = 0;
switch flag
case 'init'
    previoust = t(1);
case ''
    h = t(end) - previoust;
    previoust = t(end);
    status = (abs(h) <= hmin);
    if( y(end) < wB || y(end) > wT )
      status = 0;
    end
case 'done'
    clear previoust
end