Solving 2nd order coupled differential equations using shooting method

I’m trying to solve these two coupled 2nd order differential equations:

with the following boundary conditions:

where $r_{*}$ is the value of the $r$ at midpoint. I was trying to solve these equations using shooting method and for $l=3$ and $t=4$, I succeeded by guessing the value of $v'(frac{-l}{2})$ and $r*$:

m = 1/2 (Tanh[v[x]/(1/3)] + 1);

rv3 = ParametricNDSolve[{r[x]^4/(Guess)^2 - r[x]^2 - 
     2 r'[x] v'[x] + (r[x]^2 - m) v'[x]^2 == 0, 
   r[x]^2 - r[x]^2 v'[x]^2 - r[x] v''[x] + 2 r'[x] v'[x] == 0, 
   r[-1.5] == 10, v[-1.5] == 4, v'[-1.5] == guess}, {r, v}, {x, -1.5, 
   1.5}, {guess, Guess}, MaxSteps -> Infinity]

Manipulate[
 Plot[Evaluate[r[guess, Guess][t] /. rv3], {t, -1.5, 1.5}], {{Guess, 
   1.1039}, 1.1, 1.2}, {{guess, -9.051}, -10, -9}]
Manipulate[
 Plot[Evaluate[v[guess, Guess][t] /. rv3], {t, -1.5, 1.5}], {{Guess, 
   1.1039}, 1.1, 1.2}, {{guess, -9.051}, -10, -9}]

but for length intervals greater $l$ than 3, for example 6, I’m having problem finding the right value for guessing. Also we have symmetry along the $x$-axis at midpoint and the derivates of $r$ and $v$ with respect to $x$ are both zero $r’=v’=0$ but for lengths greater than 3, I keep getting solutions that doesn’t respect the symmetry and doesn’t have those derivatives zero at midpoint.
I’m trying to get some results similar to these:

I could replicate the $l=3$ one but for the rest, especially $l=6$ and greater I’m having problems finding the right values, because I’m getting solutions that are not correct. For example something like this for the $v-x$ plot:

where I’ve chose r[-6]==100.

Can anyone point me in the right direction so I can find the values for the equations? any help would be appreciated.