We develop numerical methods for solving wave equation with van der Pol type nonlinear boundary conditions under the framework of weak solution. Based on the wave reflection on the boundaries, we first solve the Riemann invariants by constructing two iteration mappings, and then show that the weak solution can be obtained by the integration of the Riemann invariants on the boundaries. If the compatible conditions are not satisfied or only hold with a low degree, a high order integration method is developed for the numerical solution. When the initial condition is sufficiently smooth and compatible conditions hold with a sufficient degree, we establish a sixth-order finite difference scheme, which only needs to solve a linear system at any given time instance. Numerical experiments are provided to demonstrate the effectiveness of the proposed approaches.