Using an Environmentally Friendly Renormalization Group we derive an ab initio universal scaling form for the equation of state for the O(N) model, y=f(x), that exhibits all required analyticity properties in the limits $xto 0$, $xtoinfty$ and $xto -1$. Unlike current methodologies based on a phenomenological scaling ansatz the scaling function is derived solely from the underlying Landau-Ginzburg-Wilson Hamiltonian and depends only on the three Wilson functions $gamma_lambda$, $gamma_phi$ and $gamma_{phi^2}$ which exhibit a non-trivial crossover between the Wilson-Fisher fixed point and the strong coupling fixed point associated with the Goldstone modes on the coexistence curve. We give explicit results for N=2, 3 and 4 to one-loop order and compare with known results.