This paper deals with a numerical analysis method based on a time integration technique and the H1 gradient method for solving nonparametric boundary shape optimization problems of domain boundaries in which boundary value problems of partial differential equations are defined. The H1 gradient method has been developed by applying the gradient method in a Hilbert space. In the H1 gradient method, the domain variation that minimizes the objective functional is obtained as a solution to a boundary value problem of a linear elastic continuum defined in the design domain loaded with traction in proportion to the negative shape gradients on the design boundary. The optimized shape is obtained as the solution to first-order ordinary differential equations with the solution of the H1 gradient method. The effectiveness of the proposed method is demonstrated through numerical examples in a heat conduction problem and a structural problem.