  In this article, a formulation of a point-collocation method in which the
unknown function is approximated using global expansion in tensor product
Bernstein polynomial basis is presented. Bernstein polynomials used in this
study are defined over general interval [a,b]. Method incorporates several
ideas that enable higher numerical efficiency compared to Bernstein polynomial
methods that have been previously presented. The approach is illustrated by a
solution of Poisson, Helmholtz and Biharmonic equations with Dirichlet and
Neumann type boundary conditions. Comparisons with analytical solutions are
given to demonstrate the accuracy and convergence properties of the current
procedure. The method is implemented in an open-source code, and a library for
manipulation of Bernstein polynomials bernstein-poly, developed by the authors.
