An adaptive finite element method is adopted to simulate the steady state coupled Schrödinger equations with a small parameter. We use damped Newton iteration to solve the nonlinear algebraic system. When the solution domain is elliptic, our numerical results with Dirichlet or Neumann boundary conditions are consistent with previous theoretical results. For the dumbbell and circular ring domains with Dirichlet boundary conditions, we obtain some new results that may be compared with future theoretical analysis.
