%% 5. Boundary-value problems

%%
% <latex>
% \setcounter{page}{48}
% </latex>

%%
% <latex>
% Now it is time to open the big door that will double the
% richness of our study of ODE\kern .5pt s.
% With second- and higher-order equations, since they need
% more than one boundary condition, we have the
% prospect of not just IVPs but also {\bf BVPs}, that is,
% {\bf boundary-value problems}.
% </latex>

%%
% <latex>
% Boundary-value problems are everywhere in the sciences.
% Broadly speaking,
% ODE\kern .5pt s whose independent variables represent time
% are usually IVPs modeling {\em dynamics.}
% ODE\kern .5pt s whose independent variables represent
% space are usually BVPs modeling
% {\em equilibria.}\footnote{Of course, there are
% exceptions.  A notable area where BVPs arise with a time
% variable is {\em control theory,} in which the engineer seeks to
% make a system attain certain function values at certain times.  The
% Application of Chapter~16 is in this category.}
% These differences are fundamental, and
% in this book, we will usually change the independent variable from
% $t$ to $x$ when working with a BVP.\ \ In
% particular,
% $y'$ and $y''$ will now denote
% $dy/dx$ and $d^{\kern .7pt 2}y/dx^2$ instead of $dy/dt$ and $d^{\kern .7pt 2}y/dt^2$.
% As mentioned in a footnote in Chapter 1, we
% will also plot solutions of BVPs in blue instead of green.
% </latex>

%%
% <latex>
% As our first example of a BVP, here is the solution of
% $$ y''=-y,\quad x\in [\kern .3pt 0,60\kern.3pt],
% ~y(0)=1,~y(60)=0.\eqno (5.1) $$
% The value of $y$ at $x=0$ is specified, but
% not its derivative.  A solution is found
% to match the value of $y$ specified at the other end of the interval.
% </latex>
ODEformats, L = chebop(0,60); L.op = @(x,y) diff(y,2) + y;
L.lbc = 1; L.rbc = 0; y = L\0; plot(y,CO,bvp), axis([0 60 -4 4])
title('Fig.~5.1.~~BVP (5.1) with Dirichlet BCs',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Here is the solution to another BVP, in which the condition at the
% right involves the derivative $y'$ rather than the value of $y$ itself,
% $$ y''=-y,\quad x\in [\kern .3pt 0,60\kern.3pt],~y(0)=1,~y'(60)=0.\eqno (5.2) $$
% A condition involving $y$ is called a {\bf Dirichlet} condition, and
% a condition involving $y'$ is called a {\bf Neumann} condition.  A
% condition that mixes
% both $y$ and $y'$ is called a {\bf Robin} condition.
% (To specify a Robin condition in Chebfun, such as 
% $y'(60) - y(60) = 1$, for example, we could write
% \verb|L.rbc = @(y) diff(y)-y-1|.)
% </latex>
L.rbc = 'neumann';
y = L\0; plot(y,CO,bvp), axis([0 60 -4 4])
title('Fig.~5.2.~~BVP (5.2) with a Neumann BC',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% The significance of boundary conditions becomes
% more striking if we look at equation (4.10)
% with $\omega^2 <0$, where one solution
% exponentially grows and the other exponentially decays.
% To illustrate, we revert to the Dirichlet condition of (5.1)
% and take $\varepsilon =0$, $\omega^2=-1$ in (4.10),
% $$ y''=y,\quad x\in [\kern .3pt 0,60\kern.3pt],~ y(0)=y(60)=1.\eqno (5.3) $$
% Note that the solution shown below is close to zero except near the
% endpoints, where it is far from
% zero in two regions known as {\bf boundary layers.}
% Here we begin to see strongly the difference between IVPs and BVPs.
% This curve in principle might have
% been the solution of an IVP, if the initial data had been
% chosen just right, but such a solution would hardly be likely
% to arise in practice.  On the contrary,
% from the curve it is clear that this
% solution is probably being driven by conditions imposed at the
% two ends.  We shall consider boundary layers in detail
% in Chapter~20.
% </latex>
L.op = @(x,y) diff(y,2) - y; L.lbc = 1; L.rbc = 1;
y = L\0; plot(y,CO,bvp), axis([0 60 -1.5 1.5])
title('Fig.~5.3.~~Nonoscillatory BVP (5.3)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% For the problem above, the characteristic roots (4.11)--(4.12)
% are $r_{\pm} = \pm 1$.  The solution we have plotted can
% be approximated very closely as
% $$ y(x) \approx e^{-x} + e^{x-60}. \eqno (5.4) $$
% See Exercise 5.1.
% </latex>

%%
% Suppose we are given a second-order constant-coefficient scalar linear
% homogeneous ODE BVP of the form (4.10).  To solve it, following
% Theorem 4.1, we would look for a function in the
% vector space spanned by $e^{r_+x}$ and $e^{r_-x}$, or 
% $e^{rx}$ and $xe^{rx}$ in the critically damped case, that
% satisfies the boundary conditions.
% Just as in the last chapter, a small variation of this approach
% applies to an inhomogeneous problem of the form
% $$ y'' + \varepsilon\kern .3pt y' + \omega^2\kern .3pt y = f(x), \eqno (5.5) $$
% where $f$ is a given forcing function.  If we can find
% a particular solution $y_{\rm p}^{}$ to (5.5), then the general
% solution is obtained by adding elements from this two-dimensional
% vector space.  When $f$ is simple, a particular solution
% can often be obtained by inspection or the method of undetermined
% coefficients.

%%
% <latex>
% Here is an example.  Suppose we put the
% linear function $f(x) = (x-20)/20$ on the right-hand side of (5.3),
% $$ y''=y+(x-20)/20,\quad x\in [\kern .3pt 0,60\kern.3pt],~ y(0)=y(60)=0.\eqno (5.6) $$
% A particular solution to this equation is
% $y_{\rm p}^{}(x) = -(x-20)/20$, and from
% here one can readily obtain a formula for the solution plotted below
% (Exercise 5.2).
% </latex>
L.op = @(x,y) diff(y,2) - y; L.lbc = 0; L.rbc = 0;
rhs = chebfun('(x-20)/20',[0 60]);
y = L\rhs; plot(y,CO,bvp), axis([0 60 -2 2])
title('Fig.~5.4.~~Inhomogeneous BVP (5.6)',FS,11)
%%
% <latex> \vskip 1.02em </latex>


%%
% For a nonlinear variation, let us 
% change the term $y$ in (5.6) to $y^3$,
% $$ y''=y^3+(x-20)/20,\quad x\in
% [\kern .3pt 0,60\kern.3pt],~ y(0)=y(60)=0.\eqno (5.7) $$
% The resulting curve is an interesting variation
% on Figure 5.4.  The zero of $y(x)$ near $x=20$, incidentally,
% is not exactly at $x=20$, more like $19.99999999985$.
L.op = @(x,y) diff(y,2) - y^3;
y = L\rhs; plot(y,CO,bvpnl), axis([0 60 -2 2])
title('Fig.~5.5.~~Nonlinear variation (5.7)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Returning to linear problems, let us explore some of the
% variety of boundary conditions that may arise.
% Here is an {\em advection-diffusion} equation:
% $$ 0.02y'' + y' + y = 0 , \quad x\in [\kern .3pt 0,1].
% \eqno (5.8) $$
% First let us solve it with two Dirichlet boundary conditions
% $y(0) = 0$, $y(1) = 1$.
% </latex>
L = chebop(0,1); L.op = @(x,y) 0.02*diff(y,2) + diff(y) + y;
L.lbc = 0; L.rbc = 1; y = L\0; plot(y,CO,bvp)
title('Fig.~5.6.~~Solution of (5.8) with two Dirichlet BCs',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Next we change the right boundary condition to Neumann form,
% $y'(1) = 1$.
% </latex>
L.rbc = @(u) diff(u)-1; y = L\0; plot(y,CO,bvp), ylim([-3 0])
title(['Fig.~5.7.~~Solution of (5.8) with Dirichlet ' ...
       'and Neumann BCs'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Here we couple the two endpoints with the
% conditions $y(0) = 1$, $y'(1) = -y'(0)$.
% </latex>
L.lbc = 1; L.rbc = []; L.bc = @(x,y) deriv(y,1) + deriv(y,0);
y = L\0; plot(y,CO,bvp)
title('Fig.~5.8.~~Solution of (5.8) with coupled BCs',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Not all ``boundary'' conditions must be based at
% the boundaries.  For example, we might seek a solution
% with $y(0) = 1$ and $\int_0^1 y(x) dx = 0$, so that 
% $y$ has mean zero.  A nonstandard condition
% like this is sometimes called a {\em side condition.}
% </latex>
L.bc = @(x,y) sum(y); y = L\0; plot(y,CO,bvp)
title(['Fig.~5.9.~~Solution of (5.8) with an ' ...
       'integral condition'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% So far we have considered second-order equations.
% With equations of third order
% or higher, there is the possibility of mixing boundary and interior
% conditions.  For example, here is the solution to a third-order
% problem with boundary conditions at three distinct
% points,  marked by dots in the plot for emphasis.
% $$ y''' + y = 1, \quad x\in [\kern .3pt 0,2], ~
% y(0) = 0, ~ y(1) = y(2) = 1. \eqno (5.9) $$
% </latex>
L = chebop(0,2); L.op = @(x,y) diff(y,3) + y;
L.lbc = 0; L.rbc = 1; L.bc = @(x,y) y(1)-1;
y = L\1; plot(y,CO,bvp)
hold on, plot(0:2,y(0:2),'.k',MS,14), hold off
title(['Fig.~5.10.~~Solution of 3rd-order eq.\ (5.9) ' ...
       'with interior condition'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Throughout this discussion we have avoided mentioning
% existence and uniqueness.  This is
% because these matters are not straightforward for
% BVPs, even linear ones, because of the possible presence of
% {\em eigenvalues.}   This will be the subject of the next chapter.
% For nonlinear BVPs there are further nonuniqueness
% possibilities that are both interesting and scientifically
% important, as we will explore in Chapters 16--18.  As an example of
% nonuniqueness in a nonlinear BVP, consider
% the solution that Chebfun finds to the nonlinear equation
% $$ {1\over 2} y'' + y^2 = 1, \quad x\in [-1,1], ~ y(-1) = y(1) = 0.
% \eqno (5.10) $$
% </latex>
N = chebop(-1,1);
N.op = @(x,y) .5*diff(y,2) + y^2; N.lbc = 0; N.rbc = 0; 
y = N\1; plot(y,CO,bvpnl)
title(['Fig.~5.11.~~One solution of the ' ...
       'nonlinear equation (5.10)'],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% If we ask Chebfun to look for a solution starting from
% the initial guess $y(x) = \cos(\pi x/2)$ rather than its
% default initial guess $y(x) = 0$, however, a different
% and equally valid solution is found, which we superimpose
% on the same plot.
% </latex>
N.init = chebfun('cos(pi*x/2)');
y = N\1; hold on, plot(y,'r'), hold off
title('Fig.~5.12.~~Another solution of (5.10)',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Does (5.10) have exactly these two solutions, or
% are there more?  There is no general theorem that
% gives an answer to such questions, but we shall
% make some progress on them in Chapter 16.
% </latex>

%%
% <latex>
% \begin{center}
% \hrulefill\\[1pt]
% {\sc Application: beam theory and the strength of spaghetti}\\[-3pt]
% \hrulefill
% \end{center}
% </latex>

%%
% <latex>
% In physics and engineering an idealized {\em string\/} can sustain
% a tension but has no stiffness --- no resistance
% to bending --- and this leads to second-order ODE\kern .5pt s.  This is where we would
% start, say, in studying the sounds produced by a guitar.  An idealized
% {\em beam}, on the other hand, has stiffness and resists bending, even
% though we still model it as a one-dimensional object.
% This leads to fourth-order ODE\kern .5pt s, and the same
% equations apply to the girders of the Eiffel tower as
% to a piece of dried spaghetti.
% </latex>

%%
% In the simplest setting, omitting
% dimensional constants and assuming infinitesimal
% deflections, the equation takes the form
% $$ y''''(x) = f(x) , \eqno (5.11) $$
% where $f(x)$ is the vertical force applied at position $x$ along a beam and
% $y(x)$ is the vertical deflection at that point.
% The derivatives of $y$ can be interpreted as follows:

%%
% <latex>
% \par\vskip .7em\par
% $y(x)$ and $y'(x)$: position and slope,
% \par\vskip .3em\par
% $y''(x)$: bending moment, whose square is energy density,
% \par\vskip .3em\par
% $y'''(x)$: vertical force,
% \par\vskip .3em\par
% $y''''(x)$: force density.
% \par\vskip .7em\par
% \noindent
% In a time-dependent problem, a nonzero value of $y''''(x)$ would
% be associated with vertical acceleration $-y''''(x)$, 
% but we are considering here the static problem.
% Equation (5.11) represents the condition that the force caused
% by the bending of the beam is balanced by an external force
% $f(x)$, so there is no acceleration. 
% </latex>

%%
% <latex>
% For example, suppose a beam of length $2$ and weight density $1$ is
% clamped at $x=\pm 1$,
% $$ y'''' = -1, \quad y(-1) = y'(-1) = y(1) = y'(1) = 0. \eqno (5.12) $$
% An image shows that the beam sags a little in the middle (the
% deflection is $1/24$).
% We plot $y$ with a thicker line to suggest stiffness.
% This and the next few calculations, incidentally, can be
% done on paper (Exercise 5.4), though that would change if
% complications were introduced such as variable coefficients.
% </latex>
L = chebop(-1,1); L.op = @(x,y) diff(y,4);
L.lbc = [0;0]; L.rbc = [0;0]; y = L\-1; d1 = norm(y,inf);
plot(y,CO,bvp,LW,3); axis([-1.2 1.2 -.5 .1])
title(['Fig.~5.13.~~Clamped beam (5.12) of length 2: deflection ' ...
        num2str(d1)],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% Next, instead of clamping the beam at $x=1$, let us fix the position
% there but not the slope.  In such a case one has what is called a 
% {\em natural boundary condition,}
% $y''(1) = 0$.  The deflection of the beam at the
% midpoint doubles, and the maximum deflection, which is now located
% to the right of the midpoint, increases by a factor of about 2.079.
% (The exact maximal deflection is $(39+ 55\sqrt{33}\kern .7pt )/4096.$)
% </latex>
L.rbc = @(y) [y; diff(y,2)]; y = L\-1; d2 = norm(y,inf);
plot(y,CO,bvp,LW,3); axis([-1.2 1.2 -.5 .1])
title(['Fig.~5.14.~~Beam with one freely supported ' ...
       'end: deflection ' num2str(d2)],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% Alternatively, we can let both ends be free.  Symmetry is restored,
% and the maximum deflection is in the middle again, with
% value $5/24$.
% </latex>
L.lbc = @(y) [y; diff(y,2)]; y = L\-1; d3 = norm(y,inf);
plot(y,CO,bvp,LW,3); axis([-1.2 1.2 -.5 .1])
title(['Fig.~5.15.~~Beam with two freely supported ' ...
       'ends: deflection ' num2str(d3)],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% This brings us to the spaghetti problem.  Here in Oxford, we bought a 
% pack of spaghetti at a local newsagent.
% With Scotch tape we taped 
% some pennies\footnote{Conveniently for our experiment, the British 2p coin
% weighs exactly the same as two pennies.}
% to the end of a piece, making a cantilever that can be modeled
% as a massless beam clamped at $x=0$ with a
% weight $w$ applied at $x=d$:
% $$
% y'''' = 0, \quad y(0) = y'(0) = y''(d) = 0, ~ y'''(d) = w.
% \eqno (5.13) $$
% </latex>
d = 2; L = chebop(0,d); L.op = @(x,y) diff(y,4);
L.lbc = [0; 0]; L.rbc = @(y) [diff(y,2); diff(y,3)-0.1];
y = L\0; d4 = norm(y,inf); plot(y,CO,bvp,LW,3)
penny = d+1i*(y(end)-.01)+.08*(chebfun('exp(1i*t)',[0, 2*pi])-1i);
hold on, fill(real(penny),imag(penny),[.72 .45 .2])
axis equal, xlim([-0.2 d+0.2]), hold off
title(['Fig.~5.16.~~Spaghetti cantilever (5.13): deflection ' ...
        num2str(d4)],FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% We could have measured these deflections for various cantilever
% lengths $d$ and numbers of pennies, but it seemed more fun to investigate,
% for a given number of pennies, 
% what is the length of cantilever at which the
% spaghetti breaks?
% With up to four pennies no breaking occurs --- at 26 cm, the
% piece of spaghetti is too short.
% But for $5,6,\dots, 16$ pennies we found these breaking lengths in cm:
% </latex>
d = [17.1 14.6 12.9 8.8 9.0 7.5 6.7 6.6 5.2 5.3 4.7 4.4];

%%
% <latex>
% \noindent
% Here is a plot of the data, showing how the breaking weight decreases
% as the length $d$ increases.
% </latex>
xx = linspace(4.7,17.5); plot(xx,80./xx,'-r'), hold on
plot(d,5:16,'.k',MS,12), hold off
xlabel('length of spaghetti',FS,9)
ylabel('number of pennies to break',FS,9)
text(10.3,10.4,'inverse',CO,'r',FS,10)
text(10.3,9,'linear fit',CO,'r',FS,10)
title('Fig.~5.17.~~Breaking a spaghetti cantilever',FS,11)
%%
% <latex> \vskip 1.02em </latex>

%%
% <latex>
% \noindent
% As the red line shows, the relationship appears to be
% inverse linear, and we explain this as follows.
% A piece of spaghetti breaks
% when its curvature $|y''(x)|$ exceeds a critical value for some $x$.
% From calculation or analysis of (5.13), one finds that as
% a function of $d$, $y$ scales as $wd^{\kern .7pt 3}$, where $w$ again
% is the weight.  This means that
% $y'$ scales as $wd^{\kern .7pt 2}$, and $y''$ scales as $wd$.
% If $wd$ is to remain bounded as 
% $d\to\infty$, we must have $w = O(d^{\kern .7pt -1})$.
% </latex>

%% 
% <latex>
% Returning from spaghetti to mathematics for a moment, it is interesting to
% note that the solutions of $u''''=0$ are
% precisely the polynomials of degree $\le 3$.  Thus the
% study of idealized beams with no body forces $f(x)$
% is equivalent to the study of
% cubic polynomials that satisfy various boundary conditions.
% This connection is famous in numerical analysis, where
% the classic wooden flexible rulers used for designing
% shapes in the aircraft and shipbuilding industries evolved into
% the mathematical
% {\em splines} used all across computational science and engineering.
% </latex>

%%
% <latex>
% \smallskip
% {\sc History.} Beam theory goes back to Euler and
% Jacob Bernoulli in the mid-18th century.  It was
% later that the explicit study of boundary-value problems
% became a field of mathematics in its own right with
% the investigations of {\em Sturm--Liouville problems} in the
% 1830s by Charles Sturm and Joseph Liouville.  
% French mathematics had a period of extraordinary
% productivity following the French Revolution
% and Napoleonic wars, with leading figures including Fourier
% (born 1768), Poisson (1781), Cauchy (1789),
% Sturm (1803), and Liouville (1809).
% \smallskip
% </latex>

%%
% <latex>
% {\sc Our favorite reference.}  A standard work among numerical
% analysts is Ascher, Mattheij, and Russell,
% {\em Numerical Solution of Boundary Value Problems for
% Ordinary Differential Equations,} SIAM, 1995 (first published in
% 1988).  Its appealing concrete style is well established by 
% Section I.2 on p.~7,  which begins,
% ``In this section we have collected 22 instances of BVPs which arise in
% a variety of application areas.''
% \smallskip
% </latex>

%%
% <latex>
% \begin{displaymath}
% \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl
% {\sc Summary of Chapter 5.}
% When an ODE is of order $d\ge 2$, it will normally have
% a $d$-dimensional solution space, making it possible
% for different boundary conditions to be specified at different points.
% In the simplest case $d=2$, we get a
% {\em two-point boundary value problem.} 
% Problems of elasticity and solid mechanics lead to
% fourth-order equations.
% Existence and uniqueness of solutions are
% not as straightforward for BVPs, even linear ones, as for IVPs.
% Nonlinear BVPs may have multiple solutions.
% \vspace{2pt}}}
% \end{displaymath}
% </latex>

%%
% <latex>
% \small\parindent=0pt\parskip=1pt
% {\em Exercise $5.1$.  BVPs and IVPs.}\footnote{This exercise is related to
% the numerical method for solving ODE BVPs known as {\em shooting}, and
% highlights some of the challenges of that method.  We shall say more
% about shooting in Chapter~16.}
% {\em (a)} Determine analytically the exact solution
% to (5.3).  How much does it differ from the approximation (5.4)?
% {\em (\kern .7pt b)} As remarked in the text, though a curve may be
% ``obviously'' the solution to a BVP, in principle it could equally
% have been the solution to an IVP.\ \ For this problem (5.3), what
% initial condition on $y'(0)$ would be needed to achieve
% the same solution?  {\em (c)} If that initial condition is
% multiplied by $0.99$, is the resulting value
% $y(60)$ less than or greater than Avogadro's number?
% \par
% {\em Exercise $5.2$.  Exact solution to $(5.6)$.}  {\em (a)} Derive a formula
% for the exact solution to (5.6).  {\em (\kern .7pt b)} Do the same for the case
% where the right-hand side is changed to $(x-20)^2/400$.
% \par
% {\em \underline{Exercise $5.3$}.  Sixth-order BVP.}
% The ODE
% $y^{(6)} - a y' + y = 0$, $x\in [-4,4]$, where $a\in [\kern .3pt 0,2\kern .3pt]$ is
% a parameter, has boundary conditions $y(-4) = y'(-4) = y''(-4)=0$,
% $y(4) = 1$, $y'(4)=y''(4) = 0$.  If $y(0) = 0$, what is $a$?
% \par
% {\em {Exercise $5.4$}.  Exact beam formulas.}  Find exact formulas
% for the solutions presented in Figures {\em (a)} 5.13,
% {\em (\kern .7pt b)} 5.15, and {\em (c)} 5.14.
% The reason for this ordering is that the solutions of
% Figures 5.13 and 5.15 are simpler because of symmetry.
% \par
% {\em {Exercise $5.5$}.  Exact spaghetti formula.}  Find an
% exact formula for the solution plotted in Figure 5.16.
% \par
% {\em {Exercise $5.6$}.  Cubic and quartic dependence.}
% It is readily confirmed that on
% an interval of length $d$, the deflections of Figures 5.13--5.15
% scale as $O(d^{\kern .7pt 4})$, whereas in Figure 5.16, as noted
% in the text, it is $O(d^{\kern .7pt 3})$.  
% Thus for large enough $d$ the clamped beam will sag more than
% the cantilever.  Yet the clamped beam is
% supported at {\em both\/} ends, so surely
% it should bend less.  Resolve this paradox.
% \par
% {\em \underline{Exercise $5.7$}.  Blasius equation.}
% A classic nonlinear BVP due to Blasius is
% $2\kern .3pt y''' + y\kern .3pt y'' = 0$,
% $y(0) = y'(0) = 0$, $y'(\infty) = 1$.
% Solve this with Chebfun, replacing $\infty$ by 10.  Plot the
% result and report the values of $y''(0)$ and $y(10) - 10$, both of
% which are good approximations to what one
% would get on $[\kern .3pt 0,\infty)$.
% (For more information see J. P. Boyd, ``The Blasius function
% in the complex plane,'' {\em Experimental Mathematics} 8 (1999),
% pp.~381--394.)
% \par
% {\em {Exercise $5.8$}.  Eliminating inhomogeneous boundary conditions.}
% Suppose we have a linear inhomogeneous BVP with inhomogeneous
% boundary data. To be specific, let us say it is a 
% second-order problem $Ly = f$, $y(0) = \alpha,$ 
% $y(1) = \beta$.  
% {\em (a)} Let $y_1^{}$ and $y_2^{}$ be two linearly independent
% solutions to the homogeneous problem and let $y_{\rm left}^{}$,
% $y_{\rm inner}^{}$, and $y_{\rm right}^{}$ be particular solutions
% with $f=\beta=0$, $\alpha=\beta=0$, and $f=\alpha= 0$,
% respectively.  Write expressions for a particular solution of
% the BVP and for the general solution.
% {\em (\kern .7pt b)} A different approach is as follows.  Let 
% $Y$ be {\em any\/} smooth function that satisfies $Y(0) = \alpha$ and
% $Y(1) = \beta$; $Y$ does not have to satisfy the differential equation.
% Show how the general solution of the BVP can be derived using $Y$.
% \par
% {\em \underline{Exercise $5.9$}.  Water droplet.}  The height of
% the surface of a water droplet satisfies $y'' = (y-1)(1+(y')^2)^{1.5}$
% with $y(\pm 1) = 0$.  What is the height at the midpoint $x=0$?
% \par
% {\em \underline{Exercise $5.10$}.  Nonlinear problems
% and Newton iteration.}  Chebfun solves linear BVPs by solving
% linear systems of equations obtained by discretizing in $x$, as
% described in Appendix~A.\ \ For
% nonlinear BVPs, there is an additional Newton iteration 
% involved, discussed in Appendix~A and also briefly in Chapter~16.
% {\em (a)} Rerun the problem of equation (5.3) and Figure 5.3 using
% {\tt tic} and {\tt toc} to measure how long the computation takes.
% (As with all timing experiments, do this several times to make
% sure the numbers have settled down.)  Plot the solution $y$ on
% a semilog scale and also report the length of the chebfun
% as in Exercise~1.6.
% {\em (\kern .7pt b)} Now do the same calculations again but with the right-hand
% side of (5.3) changed from $y$ to $y^{1.01}$,
% and superimpose the new solution curve
% $y$ on the same plot.  How is the computing time affected?  How is
% the length of the chebfun affected?
% {\em (c)} Show the outputs that result when {\tt L} is executed
% without a semicolon for these two computations.  Note the indication
% that in part {\em (\kern .7pt b)}, Chebfun has determined that the operator
% is nonlinear and must be treated by different numerical methods.
% </latex>
