We present the lattice simulation of the renormalization group flow
in the $3$-dimensional $O(N)$ linear sigma model.
This model possesses a nontrivial infrared fixed point, called Wilson--Fisher fixed point.
Arguing that the parameter space of running coupling constants can be spanned by
expectation values of operators evolved by the gradient flow,
we exemplify a scaling behavior analysis based on the gradient flow
in the large $N$ approximation at criticality.
Then, we work out the numerical simulation of the theory with finite $N$.
Depicting the renormalization group flow along the gradient flow,
we confirm the existence of the Wilson--Fisher fixed point non-perturbatively.