We present an algorithm for finding high order numerical approximations of minimal surfaces with a fixed boundary. The algorithm employs parametrization by high order polynomials and a linearization of the weak formulation of the Laplace-Beltrami operator to arrive at an iterative procedure to evolve from a given initial surface to the final minimal surface. For the steady state solution we measure the approximation error in a case where the exact solution is known (the catenoid). In the framework of parametric interpolation, the choice of interpolation points (mesh nodes) is directly affecting the approximation error, and we discuss how to best update the mesh on the evolutionary surface such that the parametrization remains smooth. In our test case we achieve exponential convergence in the approximation of the minimal surface as the polynomial degree increases, but the rate of convergence greatly differs with different choices of mesh update algorithms.