==> advection_LW_pbc.m <==
function [h,k,error] = advection_LW_pbc(m)
%
% Solve u_t + au_x = 0  on [ax,bx] with periodic boundary conditions,
% using the Lax-Wendroff method with m interior points.
%
% Returns k, h, and the max-norm of the error.
% This routine can be embedded in a loop on m to test the accuracy,
% perhaps with calls to error_table and/or error_loglog.

==> bvp_2.m <==
% bvp_2.m 
% second order finite difference method for the bvp
%   u''(x) = f(x),   u'(ax)=sigma,   u(bx)=beta
% Using 3-pt differences on an arbitrary nonuniform grid.
% Should be 2nd order accurate if grid points vary smoothly, but may
% degenerate to "first order" on random or nonsmooth grids. 
%
% Different BCs can be specified by changing the first and/or last rows of 
% A and F.

==> bvp_4.m <==
% bvp4.m 
% second order finite difference method for the bvp
%   u''(x) = f(x),   u'(ax)=sigma,   u(bx)=beta
% fourth order finite difference method for the bvp
%   u'' = f,   u'(ax)=sigma,   u(bx)=beta
% Using 5-pt differences on an arbitrary grid.
% Should be 4th order accurate if grid points vary smoothly.
%
% Different BCs can be specified by changing the first and/or last rows of 
% A and F.

==> bvp_spectral.m <==
% bvpspectral.m 
% pseudospectral finite difference method for the bvp
%   u''(x) = f(x),   u'(ax)=sigma,   u(bx)=beta
% The use of Chebyshev grid points gives a very accurate method.
% With a uniform grid it is quite accurate for small m, but unstable as
% m increases.
%
% Different BCs can be specified by changing the first and/or last rows of 
% A and F.

==> chap1example1.m <==
% Example 1.1 
% Table 1.1 and Figure 1.2 of  http://www.amath.washington.edu/~rjl/fdmbook

==> chap2fig2.m <==
% chap2fig2.m 
%
% second order finite difference method for the bvp
%   u'' = f, u'(ax)=sigma, u(bx)=beta
% Using 5-pt differences on an arbitrary grid.
%
% Produces Figure 2.2

==> chebpoly1.m <==
function Tnp = chebpoly1(n,z)
% evaluate first derivative of Chebyshev polynomial T_n'(z)

==> chebpoly2.m <==
function Tnpp = chebpoly2(n,z)
% evaluate second derivative of Chebyshev polynomial T_n''(z)

==> chebpoly.m <==
function Tn = chebpoly(n,z)
% evaluate Chebyshev polynomial T_n(z)

==> decay1.m <==
function decay1(tol)
% decay1.m
% Sample code for solving a system of ODEs in matlab.
% Solves  system arising from decay process  A --> B --> C.

==> decaytest.m <==
% decaytest
% test decay1 for various tolerances

==> error_loglog.m <==
function error_loglog(h,E)
% Produce log-log plot of E vs. h.
% Estimate order of accuracy by doing a linear least squares fit.

==> error_table.m <==
function error_table(h,E)
% Print out table of errors, ratios, and observed order of accuracy.

==> fdcoeffF.m <==
function c = fdcoefF(k,xbar,x)
% Compute coefficients for finite difference approximation for the
% derivative of order k at xbar based on grid values at points in x.
%
% This function returns a row vector c of dimension 1 by n, where n=length(x),
% containing coefficients to approximate u^{(k)}(xbar), 
% the k'th derivative of u evaluated at xbar,  based on n values
% of u at x(1), x(2), ... x(n).  
%
% If U is a column vector containing u(x) at these n points, then 
% c*U will give the approximation to u^{(k)}(xbar).
%
% Note for k=0 this can be used to evaluate the interpolating polynomial 
% itself.
%
% Requires length(x) > k.  
% Usually the elements x(i) are monotonically increasing
% and x(1) <= xbar <= x(n), but neither condition is required.

==> fdcoeffV.m <==
function c = fdcoeffV(k,xbar,x)
% Compute coefficients for finite difference approximation for the
% derivative of order k at xbar based on grid values at points in x.
%
% WARNING: This approach is numerically unstable for large values of n since
% the Vandermonde matrix is poorly conditioned.  Use fdcoeffF.m instead,
% which is based on Fornberg's method.
%
% This function returns a row vector c of dimension 1 by n, where n=length(x),
% containing coefficients to approximate u^{(k)}(xbar), 
% the k'th derivative of u evaluated at xbar,  based on n values
% of u at x(1), x(2), ... x(n).  
%
% If U is a column vector containing u(x) at these n points, then 
% c*U will give the approximation to u^{(k)}(xbar).
%
% Note for k=0 this can be used to evaluate the interpolating polynomial 
% itself.

==> fdstencil.m <==
function fdstencil(k,j)
% Compute stencil coefficients for finite difference approximation 
% of k'th order derivative on a uniform grid.  Print the stencil and
% dominant terms in the truncation error.
%
% j should be a vector of indices of grid points, values u(x0 + j*h)
% are used, where x0 is an arbitrary grid point and h the mesh spacing.
% This routine returns a vector c of length n=length(j) and the
% k'th derivative is approximated by
%   1/h^k * [c(1)*u(x0 + j(1)*h) + ... + c(n)*u(x0 + j(n)*h)].
% Typically j(1) <= 0 <= j(n) and the values in j are monotonically
% increasing, but neither of these conditions is required.
% The routine fdcoeffF is used to compute the coefficients.

==> heat_CN.m <==
function [h,k,error] = heat_CN(m)
% heat_CN.m
%
% Solve u_t = kappa * u_{xx} on [ax,bx] with Dirichlet boundary conditions,
% using the Crank-Nicolson method with m interior points.
%
% Returns k, h, and the max-norm of the error.
% This routine can be embedded in a loop on m to test the accuracy,
% perhaps with calls to error_table and/or error_loglog.

==> makeplotBL.m <==
% makeplotBL.m
%
% call plotBL to plot the boundary locus of the stability region S
% for a specific Linear Multistep Method.

==> makeplotS.m <==
% makeplotS.m
%
% call plotS to plot the region of absolute stability S 
% (and the Order Star if plotOS=1 below) for a specific method.   

==> odesample.m <==
function odesample(tol)
% odesample.m
% Sample code for solving a system of ODEs in matlab.
% Solves  v'' = v^2 + (v')^2 - v -1  with v(0)=1, v'(0)=0
% with true solution v(t) = cos(t).
% Rewritten as a first order system.

==> odesampletest.m <==
% odesampletest
% test odesample for various tolerances

==> plotBL.m <==
function plotBL(rho,sigma,axisbox)
% plot the boundary locus for the region of absolute stability of a Linear
% Multistep Method with characteristic polynomials rho and sigma.
%
% axisbox = axis region 
%           (if omitted, default values are used that may be a poor choice)
%
% See also makeplotBL.m, from which this routine is called for a variety of
% specific LMMs.

==> plotS.m <==
function plotS(R,plotOS,axisbox)
%
% plot S, the region of absolute stability for the rational function R(z).
%
% If plotOS=1 then also plot the Order Star (complement of the 
% region of relative stability) 
%
% Typically R(z) is an approximation to exp(z) obtained by applying a 
% finite difference method to y' = lambda y, where z = dt*lambda.
%
% axisbox = axis region 
%           (if omitted, default values are used that may be a poor choice)
%
% See also makeplotS.m, from which this routine is called for a variety of
% specific R(z) functions.

==> plotSrkc.m <==
% plotSrkc.m
%
% call plotS to plot the region of absolute stability S 
% for a Runge-Kutta-Chebyshev method of order 1 or 2.
%
% Requires
%     chebpoly.m, chebpoly1.m, chebpoly2.m
%     plotS.m from chapter7.

==> poisson.m <==
% poisson2.m  -- solve the Poisson problem u_{xx} + u_{yy} = f(x,y)
% on [a,b] x [a,b].  
% 
% The 5-point Laplacian is used at interior grid points.
% This system of equations is then solved using backslash.

==> springmass.m <==
function springmass
% springmass.m
% Oscillating spring-mass system
% Solves
%   z'' = (K/m)*(zstar - z - L) - g - (D/m) z'
% where z is the location of the mass and zstar is the location of the
% support (which may be a prescribed function of t)

==> xgrid.m <==
function x = xgrid(ax,bx,m,gridchoice)
% Specify grid points in space for solving a 2-point boundary value problem
% or time-dependent PDE in one space dimension.
%
% Grid has m interior points on the interval [ax, bx].  
% gridchoice specifies the type of grid (see below).
% m+2 grid points (including boundaries) are returned in x.