ex29ch2.m

these codes are for solving OED with matlab

function ex29ch2(fun,window,npts)
% This is a modification of dfs.m to use an output function
% instead of an event function to deal with going out the
% bottom or top of the plot window.

% 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

%>> ex29ch2('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 output function to terminate the integration on going
% out the bottom or top of the window. Also use the 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 status = outfcn(t,y,flag)
global wB wT hmin
persistent previoust
status = 0;
switch flag
case 'init'
    previoust = t(1);
case ''
    h = t(end) - previoust;
    previoust = t(end);
    status = (abs(h) <= hmin) | (y(end) < wB) | (y(end) > wT);
case 'done'
    clear previoust
end