Using MDBM to find multiple solutions of BVP problem

[seymuromarov]

I'm solving the following BVP:

using DiffEqBase, OrdinaryDiffEq
using NLsolve, Sundials
using BoundaryValueDiffEq

# ODE solver
function ode_solve(dy,y,p,t)
  #myu - gravitational parameter
  r=sqrt(y[1]^2 + y[2]^2 + y[3]^2)

  dy[1] = y[4]
  dy[2] = y[5]
  dy[3] = y[6]
  dy[4] = -myu*y[1]/r^3
  dy[5] = -myu*y[2]/r^3
  dy[6] = -myu*y[3]/r^3
end

function problem(y0,init_val,t0,t1)
    tspan = (t0,t1)
    function bc!(resid,sol,p,t)  
      resid[1] = sol[1][1]   - init_val[1]
      resid[2] = sol[1][2]   - init_val[2]
      resid[3] = sol[1][3]   - init_val[3]
      resid[4] = sol[end][1] - init_val[4]
      resid[5] = sol[end][2] - init_val[5]
      resid[6] = sol[end][3] - init_val[6]
    end

  ### Now use the BVP solver to get closer

    bvp = TwoPointBVProblem(ode_solve, bc!, y0, tspan)
    sol = solve(bvp, Shooting(Vern9()),abstol=1e-10,reltol=1e-10)

  return [sol[1][4],sol[1][5],sol[1][6],sol[end][4],sol[end][5],sol[end][6]] #kpc/Myr
end

It gives only 1 solution. However, there are at least 2 solutions.

May I implement MDBM here?