function [t, y] = w(x, I, varargin)
% W -- Compute the function w_n(t) = (t - x_0)* ... *w(t - x_n)
%
% SYNOPSIS:
%   [t, y] = w(x, I);
%   [t, y] = w(x, I, N);
%
% PARAMETERS:
%   x - Prescribed nodes x_i
%       (stored as x(1) through x(n+1) for some `n').
%   I - Interval in which the nodes are defined.  Assumed to be a
%       two-element vector, e.g. I = [-1, 1].
%
%       We also assume
%
%               I(1) <= min(x) < max(x) <= I(2).
%   N - Number of subintervals into which I is divided when sampling
%       w_n(t).
%       This parameter is OPTIONAL.  Default value is N = 100.
%
% RETURNS:
%   t - Vector of sampling points.  Equally spaced on N subintervals.
%       Thus, NUMEL(t) == N + 1.
%   y - Vector of function values for w_n.
%       y(j) = w_n(t(j))  for  j = 1:numel(t)
%
% EXAMPLE:
%   Plot w_21(t) (i.e. n = 21) on the Chebyshev nodes in [-1, 1] using
%   50 subintervals for sampling w_n(t).
%
%      % Define data
%      n = 21;
%      i = 0 : n;
%      x = cos((2*i + 1)*pi / (2*n + 2));
%
%      % Compute w_n(t)
%      [t, y] = w(x, [-1, 1], 50);
%
%      % Plot result
%      plot(t, y, 'b-', 'LineWidth', 2)
%      grid on
%
% SEE ALSO:
%   NUMEL, LINSPACE.

if nargin > 2,
   N = varargin{1};
else
   N = 100;
end

t = linspace(I(1), I(2), N + 1) .';
m = numel(x);

y = t - x(1);
for i = 2:m,
   y = y .* (t - x(i));
end

% The above loop can also be written more succinctly as the two
% vectorised statements
%
%   [T, X] = meshgrid(t, x);
%   y      = prod(T - X) .';
%
% This form is (usually) more efficient at run time than the explicit
% FOR loop, but does on the other hand require more memory.
%
% See PROD and MESHGRID for additional details.
