We present a general approach for determining the equilibrium shape of isolated, coherent, misfitting particles by minimizing the sum of elastic and interfacial energies using a synthesis of finite element and optimization techniques. The generality derives from the fact that there is no restriction on the initial or final shape, or on the elastic moduli of the particle and matrix, or on the nature of misfit. The particle shape is parameterized using a set of design variables which are the magnitudes of vectors from a reference point inside the particle to points on the particle/matrix interface. We use a sequential quadratic programming approach to carry out the optimization. Although this approach can be used to find equilibrium shapes of particles in three dimensional systems, we have presented the details of our formulation for two dimensional systems under plane strain conditions. Using our formulation, we have studied the equilibrium shapes in two dimensional systems with cubic anisotropy; the precipitate and matrix phases may have different elastic moduli, and the misfit may be dilatational or non-dilatational. The equilibrium shapes and their size dependence are analysed within the framework of symmetry-breaking shape transitions. These transitions are further characterized in terms their dependence on the cubic elastic anisotropy parameter, defined by A = 2C44/(C11 – C12), and on the modulus mismatch, defined by δ=μp/μm, where /μp and μm are the effective shear moduli of the precipitate and matrix phases, respectively. Depending on the type of misfit, the systems may be classified into the following four cases: Case A: For dilatational misfit, the equilibrium shapes in isotropic systems are circular (with an isotropic or I symmetry) at small sizes and undergo a transition at a critical size to become ellipse-like (with an orthorhombic or O symmetry). This I --O transition is continuous and is obtained only when the precipitate phase is softer than the matrix. These results are in good agreement with the analytical results of Johnson and Cahn. In cubic systems with dilatational misfit, the particles exhibit a transition from square-like shapes (with a tetragonal or T symmetry) at small sizes to rectangle-like shapes (with an O symmetry) at large sizes. This T -- O transition is continuous. It occurs even in systems with stiffer precipitates; however, it is forbidden for systems with δ >δC, where δ C represents a critical modulus mismatch. The critical size decreases with increasing cubic anisotropy (i.e., with increasing values of (A-1)/(A+1). The sides of the square-like and rectangle-like shapes are along the elastically soft directions. Case B: In these systems, the principal misfits e*xx and e*yy differ in magnitude but have the same sign. The precipitates at small sizes become elongated along the direction of lower misfit; this shape has an O symmetry. In systems with A > 1, they continue to become more elongated along the same direction, exhibiting no symmetry-breaking transition. However, in systems with A < 1, particles at large sizes are elongated along an intermediate direction between the direction of lower misfit and one of the elastically soft <11> directions; this shape has only a monoclinic or M symmetry. This O - M transition, in which the mirror symmetries normal to the x and y axes are lost, may be discontinuous or continuous. The critical size increases with δ (in the range 0.8 < δ <1.25), indicating that this transition would also be forbidden for systems with δ > δC. In systems with A < 1, the critical size decreases with increasing values of A-1/ A+1 Case C: In these systems, the principal misfits differ in both magnitude and sign, and the misfit strain tensor allows an invariant line along which the normal strain is zero. The precipitates at small sizes are elongated along the direction of lower absolute misfit, and possess an 0 symmetry. At large sizes, the mirror symmetries normal to the x and y axes are broken to yield shapes which are elongated along a direction between that of lower misfit and the invariant line. This 0 -> M transition is continuous and occurs in all the systems irrespective of the value of A The critical size increases with A and decreases with δ. Case D; The misfit in this case is a special form of that in Case C; the principal misfits have the same magnitude but opposite signs. The precipitates at small sizes have a square-like shape with its sides normal to the < 11 > axes, irrespective of the type of cubic anisotropy. At large sizes, they become rectangle-like with the long axis oriented along one of the <11> directions. Similar to Case C, this T - 0 transition is continuous and occurs in all the systems irrespective of the values of A. The critical size increases with A and decreases with δ. Thus, we have identified all the symmetry-breaking transitions in equilibrium shapes of coherent precipitates in two dimensional systems. We have identified their origin and nature, and characterized them in terms of their dependence on the anisotropy parameter and modulus mismatch.