function [M,D,K] = ccdamp(n, k)

% function [M,D,K] = ccdamp(n, k)
% 
% Spectral discretization of damped wave equation u_{tt} = u_{xx} - 2a(x) u_t
% on [0,1] with homogeneous Dirichlet boundary conditions.  The damping function
% is a(x) = 1/(x+1/k).  Castro and Cox (SIAM J Control, 2001) show that the
% spectral abscissa of the associated quadratic eigenvalue problem goes to -oo 
% as k->oo.
%
% n: discretization parameter (default: n=32).
% k: the parameter in the damping function, a(x) = 1/(x+1/k) (default: k=1e3).
% M, D, K:  Mass, damping, and stiffness matrices
%
% Uses Trefethen's "cheb.m" routine (http://www.comlab.ox.ac.uk/nick.trefethen/cheb.m )

 if nargin<1, n = 32; end
 if nargin<2, k = 1e3; end

 [Dx,x] = cheb(n);                       % Trefethen's cheb.m for [-1,1]
 Dx = Dx*2; x = (x+1)/2;                 % scale domain to [0,1]

 M = eye(n-1);                           % mass matrix
 D = 2*diag(1./(x(2:n)+1/k));            % damping matrix
 K = -Dx^2; K = K(2:n,2:n);              % stiffness matrix