This paper presents a method for optimal motion planning of continuum robots by employing Bernstein surfaces to approximate the system's dynamics and impose complex constraints, including collision avoidance. The main contribution is the approximation of infinite-dimensional continuous problems into their discrete counterparts, facilitating their solution using standard optimization solvers. This discretization leverages the unique properties of Bernstein surface, providing a framework that extends previous works which focused on ODEs approximated by Bernstein polynomials. Numerical validations are conducted through several numerical scenarios. The presented methodology offers a promising direction for solving complex optimal control problems in the realm of soft robotics.