All these functions are released under the BSD-new license.

If you find any errors or can suggest some optimisations, let me know!

Note: I will try to keep the function headers on this page in sync with the Matlab files but I can't guarantee it.

misc.zip: all the miscellaneous m-files.

expanddim.m

  function MM = expanddim(M, dim)
 
  Take an interpolation or differentiation matrix and expand it so that it
  works on multi-dimensional data sets.
 
  Written by David A.W. Barton (david.barton@bristol.ac.uk) 2008

jacobian.m

  function J = jacobian(fun, sol, [h])
 
  Construct a central difference approximation to the Jacobian of function
  fun about the argument sol.
 
  Default h = 1e-6
 
  Written by David A.W. Barton (david.barton@bristol.ac.uk) 2008

simplespectrogram.m: requires windowing routines.

  function S = simplespectrogram(data, window, overlap, shift)
 
  Calculate a short-time FFT of the data using overlapping Hamming windows.
  Data is truncated if it is not correctly related to the window size.
 
  data     Data to process (should be a 1d vector)
  window   Width of the windows to use (default length(data)/8)
  overlap  Size of the overlap between windows (default 50%)
  shift    Whether to apply fftshift to the FFT results (default true)
 
  For best performance, the window size should be a power of two. Result is
  a matrix where the columns are the frequency direction and the rows are
  the time direction.
 
  Written by David A.W. Barton (david.barton@bristol.ac.uk) 2009

integration.zip: all the numerical integration m-files. Some of these require the polynomial and miscellaneous routines as well.

odegauss.m

  function sol = odegauss(ncoll, func, trange, y0, odeopt)
 
  Integrate an ODE with an implicit Runge-Kutta collocation method at Gauss
  points (order 2, 4 or 6). The (potentially nonlinear) system is solved
  using a damped Newton iteration.
 
    ncoll     Number of collocation points (1, 2 or 3)
    func      Right-hand-side of the ODE in the form func(t, y)
    trange    Start and end times for integration. (If integration output
              at specific time points is required, use ODEGAUSSEVAL
              afterwards.)
    y0        Starting point for integration.
    odeopt    Options structure returned by ODESET. Note: only the
              InitialStep, Jacobian and Mass properties are used.
 
  Written by David A.W. Barton (david.barton@bristol.ac.uk) 2011

odegausseval.m

  function y = odegausseval(sol, x)
 
  Evaluate the solution obtained with ODEGAUSS at new time points.
 
    sol      Previously obtained solution.
    x        New time points.
 
  Written by David A.W. Barton (david.barton@bristol.ac.uk) 2011

odegauss2.m: calls ODEGAUSS with ncoll = 1 (i.e., 2nd order implicit mid-point method)

odegauss4.m: calls ODEGAUSS with ncoll = 2 (i.e., 4th order)

odegauss6.m: calls ODEGAUSS with ncoll = 3 (i.e., 6th order)

polynomials.zip: all the polynomial m-files.

bary_diff_matrix.m

  function [D M] = bary_diff_matrix(xnew, xbase, w)
 
  Calculate the both the derivative and plain Lagrange interpolation matrix
  using Barycentric formulae from Berrut and Trefethen, SIAM Review, 2004.
 
  xnew     Interpolation points
  xbase    Base points for interpolation
  w        Weights calculated for base points
 
  Written by David A.W. Barton (david.barton@bristol.ac.uk) 2008, 2011

bary_interp_matrix.m

  function M = bary_interp_matrix(xnew, xbase, w)
 
  Calculate the Lagrange interpolation matrix using Barycentric formulae
  from Berrut and Trefethen, SIAM Review, 2004.
 
  xnew     Interpolation points
  xbase    Base points for interpolation
  w        Weights calculated for base points
 
  Written by David A.W. Barton (david.barton@bristol.ac.uk) 2008, 2011

bary_weights_cheb1.m

  function w = bary_weights_cheb1(n)
 
  Calculate Barycentric weights for the n degree Chebyshev polynomial of
  the first kind.
 
  Uses explicit formulae. See Berrut and Trefethen, SIAM Review, 2004.
 
  Written by David A.W. Barton (david.barton@bristol.ac.uk) 2009

bary_weights_cheb2.m

  function w = bary_weights_cheb2(n)
 
  Calculate Barycentric weights for the n degree Chebyshev polynomial of
  the second kind.
 
  Uses explicit formulae. See Berrut and Trefethen, SIAM Review, 2004.
 
  Written by David A.W. Barton (david.barton@bristol.ac.uk) 2009

bary_weights_equispaced.m

  function w = bary_weights_equispaced(n)
 
  Calculate Barycentric weights for the n degree polynomial with base
  points equispaced.
 
  Uses explicit formulae. See Berrut and Trefethen, SIAM Review, 2004.
 
  Written by David A.W. Barton (david.barton@bristol.ac.uk) 2009

bary_weights.m

  function w = bary_weights(x)
 
  Calculate Barycentric weights for the base points x.
 
  Note: this algorithm may be numerically unstable for high degree
  polynomials (e.g. 15+). If using linear or Chebyshev spacing of base
  points then explicit formulae should then be used. See Berrut and
  Trefethen, SIAM Review, 2004.
 
  Written by David Barton (david.barton@bristol.ac.uk) 2008

poly_legendre.m

  function P = poly_legendre(n)
 
  Return the n-th Legendre polynomial.
 
  Written by David A.W. Barton (david.barton@bristol.ac.uk) 2009, 2010

poly_legendre_roots.m

  function root = poly_legendre_roots(n)
 
  Return the roots of the n-th degree Gauss-Legendre polynomial on the
  interval [-1, 1].
 
  Written David A.W. Barton (david.barton@bristol.ac.uk) 2009, 2010

quadrature.zip: all the quadrature m-files.

quad_gauss_legendre.m

  function [w, x] = quad_gauss_legendre(n)
 
  Quadrature weights for integrating f(x) over [-1, +1] at the Gauss nodes
  (i.e., the roots of the appropriate Legendre polynomial).
 
  Written by David A.W. Barton (david.barton@bristol.ac.uk) 2010, 2011

windowing.zip: all the windowing m-files.

hammingwindow.m

  function w = hammingwindow(n)
 
  Return a Hamming window of width n as given by
    w(i) = 0.54 - 0.46*cos(2*pi*i/(n - 1))
  for 0 <= i <= n - 1.
 
  Written by David A.W. Barton (david.barton@bristol.ac.uk) 2009

hannwindow.m

  function w = hannwindow(n)
 
  Return a Hann window of width n as given by
    w(i) = 0.5 - 0.5*cos(2*pi*i/(n - 1))
  for 0 <= i <= n - 1.
 
  Written by David A.W. Barton (david.barton@bristol.ac.uk) 2009