A non-perturbative Renormalization Group approach is used to calculate scaling functions for an O(4) model in d=3 dimensions in the presence of an external symmetry-breaking field. These scaling functions are important for the analysis of critical behavior in the O(4) universality class. For example, the finite-temperature phase transition in QCD with two flavors is expected to fall into this class. Critical exponents are calculated in local potential approximation. Parameterizations of the scaling functions for the order parameter and for the longitudinal susceptibility are given. Relations from universal scaling arguments between these scaling functions are investigated and confirmed. The expected asymptotic behavior of the scaling functions predicted by Griffiths is observed. Corrections to the scaling behavior at large values of the external field are studied qualitatively. These scaling corrections can become large, which might have implications for the scaling analysis of lattice QCD results.