Numerical solutions of surface and internal solitary waves are obtained through a new method, where advection equations on physical quantities including surface/interface displacements and velocity potential are solved to find convergent solutions by applying the Newton-Raphson method. The nonlinear wave equations derived using a variational principle are adopted as the fundamental equations in the present study. Surface and internal solitary waves obtained through the proposed method are compared with the corresponding theoretical solutions, as well as numerical solutions of the Euler equations, to verify the accuracy of solutions through the proposed method especially for internal solitary waves of large amplitude progressing at a large celerity with a flattened wave profile. Properties of surface and internal solitary waves are discussed considering vertical distribution of horizontal and vertical velocity, as well as kinetic and potential energy. Numerical simulation using a time-dependent model has also been performed to represent propagation of surface and internal solitary waves with permanent waveforms.