Stiff second order ODE

The equation is not stiff, despite the claim by Mathematica. It can be solved by using the "Shooting" Method.

NDSolve[{f''[x] == 2 (f[x]^3 - f[x]), f[-3] == 1, f[3] == -1}, f, {x, -3, 3}, 
 Method -> {"Shooting", "StartingInitialConditions" -> {f[-3] == 1, f'[-3] == 0}}][[1, 1]];
Plot[f[x] /. %, {x, -3, 3}, AxesLabel -> {f, x}]

This looks like a plot of -Tan[x] but is not precisely that because of the boundary conditions.

Addendum - Sample Solution for Question in Comment

satoru describes a much more difficult problem in a comment below. A sample solution to it is

xm = 10;
NDSolve[{f''[x] == f[x]^3 - f[x] + α f[x] (1 - Tanh[x])^2, 
    f[-xm] == Sqrt[1 - 2 α], f[xm] == -1} /. α -> 10^-6, f, {x, -xm, xm}, 
    WorkingPrecision -> 30, Method -> {"Shooting", 
    "StartingInitialConditions" -> {f[xm] == -1, f'[xm] == 0}}][[1, 1]];
Plot[f[x] /. %, {x, -xm, xm}, AxesLabel -> {f, x}]

Note how the Tanh-like solution is shifted to the left even for tiny α. Solutions for larger a probably can be obtained by centering the integration range on the value of x for which f[x] == Sqrt[1 - 2 α] - 1, and obtaining an estimate for that quantity by extrapolating results from smaller values of α. See, for instance, the automated process in my answer here.

ClearAll["Global`*"];
Remove["Global`*"];

sol = NDSolve[{f''[x] == f[x]^3 - f[x], f[-3] == 1, f[3] == -1}, f[x], {x, -3, 3}, 
Method -> {"Shooting", "StartingInitialConditions" -> {f'[0] == 0, f[0] == 1}}];

UPDATE:

Maples symbolic solution(The symbolic solution is quite long and complex) ,converted to numeric form.Maple gives 5 solutions other than Mathematica.

Maples JacobiSN[z,k] is equal to Mathematica JacobiSN[z,k^2].

With 20 digits precision.

maple1 = (-1.5279997865913399554 - 0.58636062345487262098*I)*
JacobiSN[(0.94894046904339260202 - 
  0.94416766565747605255*I)*(0.70710678118654752440*x + 
  1.6985155326343786960 - 
  0.20793350722877537103*I), (-0.50021641106619559220 - 
  1.1156113783217666782*I)^2]; 
maple2 = (-0.57150532368340334446 - 
0.87200104339127158745*I)*
JacobiSN[(1.5911818830635583675 - 
  0.31319690342130914617*I)*(0.70710678118654752440*x + 
  1.7471798211151876014 + 
  0.44795429731942814652*I), (-0.24192870056640322887 - 
  0.59564049424232431474*I)^2];
maple3 = (-1.5279997865913399554 - 
0.58636062345487262098*I)*
JacobiSN[(0.94894046904339260202 - 
  0.94416766565747605255*I)*(0.70710678118654752440*x + 
  1.6985155326343786960 - 
  0.20793350722877537103*I), (-0.50021641106619559220 - 
  1.1156113783217666782*I)^2]; 
maple4 = (-1.5279997865913399554 - 
0.58636062345487262098*I)*
JacobiSN[(0.94894046904339260202 - 
  0.94416766565747605255*I)*(0.70710678118654752440*x + 
  1.6985155326343786960 - 
  0.20793350722877537103*I), (-0.50021641106619559220 - 
  1.1156113783217666782*I)^2]; 
maple5 = (-1.5279997865913399554 + 
0.58636062345487262098*I)*
JacobiSN[(0.94894046904339260202 + 
  0.94416766565747605255*I)*(0.70710678118654752440*x + 
  1.6985155326343786960 + 
  0.20793350722877537103*I), (-0.50021641106619559220 + 
  1.1156113783217666782*I)^2];

.

Boundary conditions check:

  Re[{maple1, maple2, maple3, maple4, maple5}] /. x -> -3 // N
  (*{1., 1., 1., 1., 1.}*)
  Re[{maple1, maple2, maple3, maple4, maple5}] /. x -> 3 // N
  (*{-1., -1., -1., -1., -1.}*)

.

  Plot[{Re[maple1], Re[maple2], Re[maple2], Re[maple2], Re[maple2], 
  f[x] /. sol}, {x, -3, 3}, 
  PlotLegends -> {"maple1", "maple2", "maple3", "maple4", "maple5", 
  "NDSOLVE"}, 
  PlotStyle -> {Red, {Green, Dashing[{0.2, 0.05}], 
  Thickness[0.01]}, {Blue, Thickness[0.01], 
  Dashing[{0.3, 0.1}]}, {Black, Thickness[0.01], 
  Dashing[{0.1, 0.1}]}, Yellow, Brown}, AxesLabel -> {x, f[x]}]

In versiion 12.0 you can use the nonlinear FEM solver:

sol = NDSolveValue[{f''[x] == f[x]^3 - f[x], f[-3] == 1, f[3] == -1}, 
  f[x], {x, -3, 3}, Method -> "FiniteElement"]
Plot[sol, {x, -3, 3}, AxesLabel -> {f, x}]