function [y,xmat] = lincmb(t,c,xfun,m,varargin)
%
% Usage:  lincmb(t,c,xfun,m)
%     or  lincmb(t,c,xfun,m,s)
%
%     where if the optional argument s is
%     omitted, it is taken to be 1.
%
% Compute
%
%   y(t_i) = \sum_j c_j xfun_j(t_i)
%
%   and
%
%   xmat = matrix with ij element xfun_j(t_i)
%
%   where xfun_j(t_i) = xfun( s_j(t_i-m_j) )
%
% If you want to plot the funtions 
% xfun(s_j(t-m_j)) themselves instead of
% their linear combination, you can
%
%      plot(t,xmat)
%
% t = ARRAY of times at which the
%     linear combination will be evaluated.
% c = ROW vector of coefficients of the
%     linear combination.
% xfun = either a STRING containing the name of
%        the underlying function to be used or
%        the handle of an anonymous function.
%
% m = ROW vector of shifts.
% s = OPTIONAL ROW vector of scale factors.
%     WARNING: If length(s)>1, then s and
%     m are assumed to have the same length;
%     otherwise an ERROR may result.
%
% y = vector output (same dimension as t)
% xmat = matrix (see below).
%
% The key idea is to observe that
% the formula
%
% y(t_i) = \sum_{j=1}^n c_j xfun(t_i-m_j)
%
% can be viewed as matrix-vector
% multiplication if we think of c and
% y as column vectors.
%
sizt = size(t); % Make t a column vector
lt = prod(sizt);
tt = reshape(t,lt,1);
lm = length(m);

if nargin==4
  tmat = repmat(tt,1,lm);
  mmat = repmat(m,lt,1);
else
  s = varargin{1};
  ls = length(s);
  if ls==1
    tmat = repmat(s(1)*tt,1,lm);
  else
    tmat = tt*s;
  end
  mmat = repmat(s.*m,lt,1);
end
xmat = feval(xfun,tmat-mmat);
yloc = xmat*(c.');
y = reshape(yloc,sizt);
