My question: is there a way to determine the "stable" region in each Lagrange location? Particularly in the L4 and L5 regions.
tl;dr: Yes but it's usually a fancy version of "trial and error".
We are sometimes wrongly told that this abominable zero-velocity potential diagram show areas that are "bound" or "unbound" to Lagrange points, but it does not and that's obvious, because while the surface dips near L1, L2 and L3 it is a maximum near L4 and L5 and yet we know that for some mass ratios these regions can collect and accumulate trojan asteroids.
Source
The problem is that these are 2D plots of the pseudo-potential at one specific velocity. The pseudo-potential is given by
$$C = \omega^2 (x^2+y^2) + 2\left(\frac{\mu_1}{r_1} + \frac{\mu_2}{r_2} \right) - (v_x^2 + v_y^2 + v_z^2)$$ and these zero velocity surfaces which are contours when drawn in a 2D slice don't depict orbits at all. If you put a particle on any contour other than $C=0$ it will accelerate rather than follow one of these lines.
So instead what people do is to run numerical simulations with a whole range of starting position and velocity vectors and find out by a fancy version of "trial and error" which particles will end up back at the same position where they started with the same velocity vector they started with. Every time they find a closed, periodic orbit it helps map out the "stable region".
One example of a method to find regions of 4D or 6D phase space ($x, y, v_x, v_y $ or $x, y, z, v_x, v_y, v_z$) that have closed periodic orbits is the shooting method which can be made faster by use of the transition matrix.
In Figure 1 of F. Marzari et al. Origin and Evolution of Trojan Asteroids you can see how tempting it is to use a zero velocity contour to indicate where the region of stable orbits lay with out actually saying it. Instead the the words "are closely related to" are used.
Fig. 1. The location of the five Lagrange equilibrium points in the circular-restricted three-body problem. The primary and secondary masses (the Sun and planet in our examples) are denoted by the large and small filled circles respectively. The selected zero velocity curves (see text) are closely related to the types of orbits that can occur in the system. The letters P, H, and T denote passing, horseshoe, and tadpole orbits respectively. Note that the two masses and each of the L4 and L5 points form an equilateral triangle.
If you want to go beyond just a circular restricted three-body model and to the real solar system, numerical simulations become much more complex because there is long term evolution of these orbits due to all of the devitaions of reality from the simple CR3BP model. See just for an example The resonant structure of Jupiter's Trojan asteroids – I. Long-term stability and diffusion
