We use bivariate splines to solve Helmholtz equation with large wave number, e.g. wave number k=1000 or larger while the size of underlying triangulation is reasonable. It is known that for linear finite elements, the wave number kand the size h of triangulation satisfies a stubborn relation k2h<1. For high order finite elements, so-called hp version of internal penalty discontinuous Galerkin (IPDG) methods, the relation is k3h2=O(p2) for polynomial degree p. For continuous internal penalty finite element method (CIP-FEM), the relation is O(k2p+1h2p)=O(k1+1/(2p)hp). That is, when k=1000, h has to be extremely small, i.e. the size of underlying triangulation has to be extremely small which leads to the infeasibility of numerical computation. In this talk, we demonstrate that the bivariate spline method enables us to have very accurate solution for kh ≥ p/2 when we use p ≥ 10. No pollution phenomenon is observed in our experiments as long as kh=p/2. We then extend our spline method to numerical solution of a Helmholtz problem over inhomogeneous media and a wave equation with periodic source arising from Maxwell's equations in potential formulation. In addition, I will explain the existence, uniqueness and stability of our spline solutions as well as their approximation.