function omega = bidim ( n, alpha, beta, d )

%*****************************************************************************80
%
%% bidim() solves a bidimensional Vandermonde system.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 November 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gustavo Galimberti, Victor Pereyra,
%    Numerical Differentiation and the Solution of Multidimensional Vandermonde
%    Systems,
%    Mathematics of Computation,
%    Volume 24, Number 110, April 1970, pages 357-364.
%
%  Input:
%
%    integer N, the order of the data vectors.
%
%    real ALPHA(N), BETA(N), the values that define A.
%
%    real D(((N+1)*N)/2), the right hand side.
%
%  Output:
%
%    real OMEGA(((N+1)*N)/2), the solution.
%
  nn = ( ( n + 1 ) * n ) / 2;
%
%  FOR NOW, SET P = NN
%
  p = n;

  c = zeros ( p, p );
  omega = zeros ( nn, 1 );

  for k = 1 : p
%
%  i) Solve a p+1-k by p+1-k Vandermonde system for C(K,1:P+1-K);
%
    if ( k == p )
      c(k,1) = d(p);
    else
      c(k,1:p+1-k) = pvand ( p+1-k, beta(1:p+1-k), d(1:p+1-k) )
    end
%
%  ii) Obtain OMEGA(1:K).
%
    if ( k == 1 )
      koff = 0;
      omega(1) = c(1,p)
    else
      koff = ( k * ( k - 1 ) ) / 2;
      omega(koff+1:koff+k) = pvand ( k, alpha(1:k), c(1:k,p+1-k) )
    end

    if ( k == p )
      break
    end

    pause
%
%  iii) Compute C(K+1:P,P+1-K) from OMEGA(1:K).
%
    for nu = k + 1 : p
      c(nu,p+1-k) = 0.0;
      for i = 1 : k
        c(nu,p+1-k) = c(nu,p+1-k) + alpha(i) ^ ( nu - 1 ) * omega(koff+i);
      end
    end
    c
    pause
%
%  iv) Update D:
%
    munu = ( ( k + 1 ) * k ) / 2;
    for nu = k + 1 : p
      for mu = 1 : p + 1 - nu
        munu = munu + 1;
        d(munu) = d(munu) - beta(p+1-k) ^ ( mu - 1 ) * c(nu,p+1-k);
      end
    end
    d 
    pause

  end

  return
end

