  The Wilsonian renormalization group (RG) method is applied to finite
temperature systems for the study of non-perturbative methods in the field
theory. We choose the O(N) linear sigma model as the first step. Under the
local potential approximation, we solve the Wilsonian RG equation as a
non-linear partial differential equation numerically. The evolution of the
domain is taken into account using the naive cut and extrapolation procedure.
Our procedure is shown to yield the correct solution obtained by the auxiliary
field method in the large $N$ limit. To introduce thermal effects, we consider
two schemes. One in which the sum of the Matsubara frequencies are taken before
the scale is introduced is found to give more physical results. We observe a
second order phase transition in both the schemes. The critical exponents are
calculated and are shown to agree with the results from lattice calculations.
