There are numerous MATLAB resources for spectral and pseudospectral methods.

(1) Gautschi, W. Algorithm 726: ORTHPOL—A Package of Routines for Generating Orthogonal Polynomials and Gauss-Type Quadrature Rules, ACM Trans. Math. Software 20, 21-62 (1994). The source codes are here.

(2) MATLAB Differention Suite by Weideman and Reddy ACM Trans. Math. Sci. 26, 465-519 (2000). The source codes are here.

(3) Trefethen, L. N. Spectral Methods in Matlab, SIAM (2000).

(4) Gautschi, W. Orthogonal Polynomials (in MATLAB), J. Comput. Appl. Math 178, 215-234 (2005).

(5) CHEBFUN suite of Matlab codes from Oxford University at chebfun.org

**MATLAB codes that accompany Spectral Methods in Chemistry and Physics.**

Chapter 3 Chapter 4 Chapter 5 Chapter 6 (Click on a chapter to be redirected.)

Chapter 2 – Click on a file and save, changing the file extension from *.txt to *.m.

ab_bimode | Calculation of recurrence coefficients, β_{n}, for full range bimodal polynomials written to bimodal.dat. Coefficients for ε = 0.02 and 0.005 are in bimodal_eps02 and bimodal_eps005 and shown graphically in beta_n_eps02 and beta_n_eps005. |

bimode_polyn_with_pts_eps02 | Calculation of full range bimodal polynomial of order n for ε = 0.02 with recurrence coefficients in bimodal_eps02.dat. sample output bimodal_polyn_with_pts_eps02. Quadrature points are shown by the solid circles. Maple code bimode-beta_n-eps02-Maple |

bimode_polyn_with_pts_eps005 | Calculation of full range bimodal polynomial of order n for ε = 0.005 with recurrence coefficients in bimodal_eps005.dat. sample output bimodal_polyn_with_pts_eps005. Quadrature points are shown by the solid circles. |

LnSqwtspts | MultiExp quadrature by Gill and Chien (2003). Sample output for N = 10, LnSqwtspts10 – see also Gill web site |

ab_maxwell_p0 | Recurrence coefficients, α_{n} and β_{n}, for Maxwell polynomials, p=0; output data abmaxp0 – more precise data from Gautschi abhrhermite |

ab_maxwell_p2 | Recurrence coefficients, α_{n} and β_{n}, for Maxwell polynomials, p=2; output data abmaxp2 |

maxwell_p0_with_pts | Maxwell polynomial for p = 0 and N = 10 in the graph Maxwell_p0_with_pts compared with Hermite (N = 20) in Fig. 2.4 for x positive. |

Rys_pw | Recurrence coefficients, α_{n} and β_{n}, for Rys polynomials, xmin = -1 for full range, J_n(x), and xmin = 0 for half range, R_n(x). For N = 13, xmin = -1, recurrence coefficients Rys_fr_ab1 and quadrature points Rys_fr_pw1 are to be compared with Sagar and Smith. For N = 40, xmin = -1, recurrence coefficients Rys_fr_ab2 and quadrature points Rys_fr_pw2 are to be compared with Schneider and Nygaard. |

Rys_Jn_with_pts | Full range, J_n(x), n = 20, c=2, polynomials as shown in Fig. 2.9 Jn_with_pts. Polynomials are not normalized to unity. |

Rys_Rn_with_pts | Half range, R_n(x), n = 10, c = 2, polynomials as shown in Fig. 2.9 Rn_with_pts. Polynomials are not normalized to unity. |

ab_GK_rad_transfer | Recurrence coefficients, α_{n} and β_{n}, for w(x)=exp(-3/2x) introduced by Gander and Karp for radiative transferproblems. Output in ab_GK is to be compared with results with Maple in Table 5 in Gander and Karp |

Legptswts | Calculation of Legendre quadrature weights and points. Sample output for N = 10, Legptswts10 – see also the Keisan online calculator |

Hermptswts | Calculation of Hermite quadrature weights and points. Sample output for N = 10, Hermptswts10 – see also the Keisan online calculator |

Lagptswts | Calculation of Laguerre quadrature weights and points. Sample output for N = 10, Lagptswts10 – see also the Keisan online calculator |

lagpoly_with_pts | Graph for the Associated Laguerre polynomial with quadrature points – click to download sample output |

hermite_with_pts | Graph for the Hermite polynomials with quadrature points – click to download sample output |

Legendre_with_pts | Graph for the Legendre polynomials with quadrature points – click to download sample output |

Chapter 3: – Click on a file and save, changing the file extension from *.txt to *.m.

legendre_norm_linear | Quadrature evaluation of the normalization of two Maxwell polynomials, n=4 and n=6. Convergence with the number of Legendre quadrature points, N, and a linear map to [-1,1] as shown in Fig. 3.3. |

legendre_norm_expmap | Quadrature evaluation of the normalization of two Maxwell polynomials, n=4 and n=6. Convergence with the number of Legendre quadrature points, N, and an exponential map to [-1,1] as shown in Fig. 3.3. |

legendre_norm_Becke | Quadrature evaluation of the normalization of two Maxwell polynomials, n=4 and n=6. Convergence with the number of Legendre quadrature points, N, and the Boyd-Becke map to [-1,1] as shown in Fig. 3.3. |

legendre_norm_MK | Quadrature evaluation of the normalization of two Maxwell polynomials, n=4 and n=6. Convergence with the number of Legendre quadrature points, N, and the Mura-Knowles map to [-1,1] as shown in Fig. 3.3. |

radint_maxp2_vsN_a | Convergence of the radial integral, Eq. (3.10), for 2 Gaussians, Eq. (3.14), versus N for several scale factors s as shown in Fig. 3.6(A). The linear fit parameters A(s) and B(s) are printed to the screen and are shown in Table 3.2. |

radint_maxp2_vsN_b | Convergence of the radial integral, Eq. (3.10), for 3 Gaussians, Eq. (3.15), versus N for several scale factors s as shown in Fig. 3.6(B). The linear fit parameters A(s) and B(s) are printed to the screen and are shown in Table 3.2. |

lindh_integ_max2 | Integration of a Gaussian, Eq. (3.11), with Maxwell (p = 2) quadratures; with results in Fig.3.4. Data file with recurrence coefficients for Maxwell quadratures is in abmaxp2 |

lindh_integ_max2B | Integration of a Gaussian, Eq. (3.11), with Maxwell (p = 2) quadratures versus the number of quadrature points, N for different values. |

sigma_O_O2 | Energy dependence of the reactive cross section, Eq. (3.25), for collisional dissociation of O-O2 reactive collisions in comparison with the line-of-centers cross section, Eq. (3.23), as shown in Fig. 3.7. |

integrand_O_O2 | Variation of the integrand for the calculation of the equilibrium O-O2 dissociative rate coefficient as shown in Fig. 3.7 for three temperatures. |

sr_rate_O_O2 | Convergence of the Simpson’s rule integration of the O-O2 dissociative rate coefficient as shown in Table 3.4. |

lag_rate_O_O2 | Convergence of the equilibrium O-O2 dissociative rate coefficient with Laguerre quadratures as shown in Table 3.4. |

lag_scaled_rate_O_O2 | Convergence of the equilibrium O-O2 dissociative rate coefficient with Laguerre quadratures as shown in Table 3.4 versus the scaling parameter, s. |

nurx_tdep_speed | This code evaluates the thermal average, Eq. (3.22), of the nuclear cross section, Eq. (3.33). The convergence of the nuclear fusion reaction rate evaluated with a Maxwell quadrature based on the weight function w(x) = xexp(-x*x). Inp ut parameters are; tk = temperature in 10**9 K, sc = quadrature scale factor and irx identifies a reaction from a choice of four. The results are in Table 3.7 |

Chapter 4: – Click on a file and save, changing the file extension from *.txt to *.m.

sin_expand_hermite | Expansion of sin(x) in Hermite polynomials; Fig. 4.4. The MATLAB code is run three times to produce this figure. Some changes by hand are required between each run to format the graph as noted. |

plot_kappa | Plot_kappa plots the Kappa distribution for 5 kappa values in comparison with the Maxwellian |

expand_kappa_lag | Expansion of the Kappa distribution in Laguerre polynomials demonstrating the nonconvergence; Fig. 4.8 for kappa = 30 with s = 1 and nmax = 25. For kappa = 50, the plotting format has to be changed. |

wave_packet_cosines | Plot the graph for the sum of cosines, Eq. (4. 105) as shown in Fig. 4.12 |

wave_packet | A model for a quantum mechanical wavepacket as shown in Fig. 4.13 |

fid | A model free induction decay curve and the Fourier transform are shown in Fig.4.14. The frequency spectrum of the signal is recovered with the FFT. |

fs_step | Fourier sine series of a step function shown in Fig. 4.15 as an introduction to Gibbs phenomenon |

runge_uniform | Demonstration of the Runge phenomenon with a uniform interpolation grid for the function, Eq. (4.151) as shown in Fig. 4.20 (upper graphs). |

runge_cheb | Interpolation of the function, Eq. (4.151) with Chebyshev quadrature points as shown in Fig. 4.20 (lower graph) |

Chapter 5: – Click on a file and save, changing the file extension from *.txt to *.m.

noneq_rx | Nonequilibrium effects for an elementary chemical reaction as shown in Fig. 5.4(A) |

noneq_spectral | Spectral convergence of the solution of the Boltzmann equation for an elementary chemical reaction as summarized in Fig. 5.4(C) |

linbe_gr_eig | Eigenvalues of the linearized Boltzmann collision operator, Eq. (5.65). Spectral convergence of the “spectral gap”, λ_{2}.Sonine-Laguerre basis functions are used. |

linbe_eig | Eigenvalues of the linearized Boltzmann collision operator with the pseudospectral method based on the Maxwell polynomial. Maxwell polynomial quadratures are used, Eqs. (5.99) and (5.100). This requires the recurrence coefficients up to N = 200 in ab_maxp2_200 as a *.dat file.(This code could be further vectorized.) |

linbe_eig_gamma1 | Eigenvalues of the linear Boltzmann collision operator, Eqs. (5.108) and (5.110), for mass ratio =1. This requires the recurrence coefficients up to n = 200 in ab_maxp2_200 as a *.dat file. (This code could be further vectorized.) |

rt_chandra | Eigenvalues of the radiative transfer equation, Eq. (5.19), and the extrapolation length. Eq. (5.25) |

Chapter 6: – Click on a file and save, changing the file extension from *.txt to *.m.

Harm_Osc_Fourier | Convergence of the eigenvalues for the quantum harmonic oscillator with Fourier basis functions as shown in Fig. 6.14(A). |

Harm_Osc_Herm1 | Eigenvalues of the quantum harmonic oscillator with Hermite basis functions with the representation as in Eq. (6.175) as in Baye and Heenen (1986). The results obtained with this Matlab code shows “ghost” levels as in Table 6.6. |

Harm_Osc_Herm2 | Exact calculation of the eigenvalues of the quantum harmonic oscillator with Hermite basis functions with the representation as in Eq. (6.178). The derivative operator is constructed as in Eq. (3.138) without explicit reference to the harmonic potential. The eigenvalues are exact as are the eigenfunctions shown in Fig. (6.15) calculated with this code. |

SE_electron_potl | The convergence of the eigenvalues for the Schroedinger equation with the potential, Eq. (6.180), based on the nonclassical polynomials orthogonal with respect to the weight function, Eq. (6.181). The recurrence coefficients for this weight function are in the file ab_se_electron. This Schroedinger equation is isospectral with the Fokker-Planck equation for electron relaxation for a hard sphere cross section, namely Eq. (6.63), with sigma(x) = 1 |