We study the curvaton scenario using gauge-invariant second order perturbation theory and solving the governing equations numerically. Focusing on large scales we calculate the non-linearity parameter fNL in the two-fluid curvaton model and compare our results with previous analytical studies employing the sudden decay approximation. We find good agreement of the two approaches for large curvaton energy densities at curvaton decay, Ωσdec, but significant differences of up to 10% for small Ωσdec.