How can I shade specific regions between two curves in a Plot?

With ParametricPlot, you can use Mesh and its friends. The MeshShading is a matrix of colors that maps onto the grid created by the mesh. Here it's 3 x 4 color matrix. Apply ParametricPlot to {s13, f (1 - t) + g t} which interpolates between the two curves.

With[{f = ((0.707106 Sqrt[1 - 1.5 s13^2] + 0.5 s13)^2)/(1 - s13^2),
      g = ((0.707106 Sqrt[1 - 1.5 s13^2] - 0.5 s13)^2)/(1 - s13^2)}, 
 Show[
  ParametricPlot[
   {s13, f (1 - t) + g t}, {s13, 0.0, 0.23}, {t, 0, 1},
   PlotRange -> {{0.0, 0.23}, {0.35, 0.67}},
   MeshFunctions -> {#1 &, #2 &},  (* x, y *)
   Mesh -> {{0.04, 0.1, 0.15}, {0.46, 0.54}}, (* x, y coords *)
   MeshStyle -> None,
   MeshShading -> {
     {Directive[Opacity[.19], Green], Directive[Opacity[.19], Green],
       Green, Directive[Opacity[.19], Green]},
     {Directive[Opacity[.19], Green], Yellow, Brown, Yellow},
     {Directive[Opacity[.19], Green], Directive[Opacity[.19], Green],
       Green, Directive[Opacity[.19], Green]}},
   Frame -> True, FrameLabel -> {sin13, sin223},
   AspectRatio -> 1/GoldenRatio],
  Plot[{f, g}, {s13, 0.0, 0.23}, 
   PlotStyle -> {Directive[Gray, Thick], Directive[Brown, Thick]}]
  ]
 ]

I had to add the boundary lines by hand using Plot, since there can be only a single BoundaryStyle.

(Edit notice: Originally, I goofed and made the boundary line straight lines.)

Perhaps something like:

Plot[{x, -x, 
       UnitStep[x - 2.5],   -UnitStep[x - 2.5], 
      x UnitBox[x - 2.5], - x UnitBox[x - 2.5]}, 
      {x, 0, 4}, 
      PlotStyle -> Join[Blue {1, 1}, ConstantArray[None, 4]],
      Filling -> {1 -> {2}, 3 -> {4}, 5 -> {6}}]