The Beyond-Linear-Use-of-Equation-Superposition (BLUES) function method is extended to coupled nonlinear ordinary differential equations and applied to the epidemiological SIRS model with vaccination. Accurate analytic approximations are obtained for the time evolution of the susceptible and infected population fractions. The results are compared with those obtained with alternative methods, notably Adomian decomposition, variational iteration and homotopy perturbation. In contrast with these methods, the BLUES iteration converges rapidly, globally, and captures the exact asymptotic behavior for long times. The time of the infection peak is calculated using the BLUES approximants and the results are compared with numerical solutions, which indicate that the method is able to generate useful analytic expressions that coincide with the (numerically) exact ones already for a small number of iterations.