Am I seeing a particle orbiting a Morris-Thorne Wormhole?

I tried to calculate the time-like geodesics of the Morris-Thorne Wormhole $[1]$, $[2]$, $[3]$ for redshift function $Phi(r) = 0$ and $b(r) = sqrt{r_{0} r}$. But I don’t know for sure if all the plots that I made are correct.

So, the metric is:

$$ds^2 = -e^{2Phi(r)}dt^2+Biggr{ 1-frac{b(r)}{r}Biggr}^{-1}dr^2+r^2[dtheta^2+sin^2(theta)dphi^2]tag{1}$$

And the plots are:

The various spheres are the 2-spheres due to spherical symmetry of the metric, and the central sphere is the sphere of the radius $r_{0}$, therefore is the throat.

My doubt is:

Am I seeing a particle orbiting the throat of a Morris-Thorne Wormhole? In other words is it correct to say that the particle approches from infinity, orbits the throat and returns to the same universe?

$$ * * * $$

$[1]$ LOBO.F.S.N; Exotic solutions in General Relativity,arxiv:0710.4474v1, 2007

$[2]$ LOBO.F.S.N; Wormholes, Warp Drives and Energy Conditions,Springer, Vol.189, 2017

$[3]$ VISSER.M; Lorentzian Wormholes, Springer, AIP Press, 1995

Appendix: Mathematica code:

I posted this question here as well, due to the mathamtica code.

n = Length[coords];
a = 0;
r0;
[CapitalPhi][r] = 0;
B[r] = (Sqrt[r*r0]);

tt = -(Exp[
    2*[CapitalPhi][
      r]   ]); rr = (1)/(1 - ((B[r])/(r))); [Theta][Theta] = 
 r^2; [CurlyPhi][CurlyPhi] = 
 r^2*(Sin[[Theta]]*Sin[[Theta]]); t[CurlyPhi] = 0; metric = {{tt, 
   0, 0, 0}, {0, rr, 0, 0}, {0, 0, [Theta][Theta], 0}, {0, 0, 
   0, [CurlyPhi][CurlyPhi]}}; metric // MatrixForm
inversemetric = Simplify[Inverse[metric]]; inversemetric // MatrixForm
christoffel := 
 christoffel = 
  Simplify[Table[
    1/2*Sum[inversemetric[[i, 
        s]]*(D[metric[[s, j]], coords[[k]]] + 
         D[metric[[s, k]], coords[[j]]] [Minus] 
          D[metric[[j, k]], coords[[s]]]), {s, 1, n}], {i, 1, n}, {j, 
     1, n}, {k, 1, n}]]
listchristoffel := 
 Table[If[UnsameQ[christoffel[[i, j, k]], 
     0], {ToString[[CapitalGamma][i, j, k]], 
     christoffel[[i, j, k]]}], {i, 1, n}, {j, 1, n}, {k, 1, 
    j}] TableForm[
   Partition[DeleteCases[Flatten[listchristoffel], Null], 2], 
   TableSpacing -> {2, 2}]
geodesic := 
 geodesic = 
  Simplify[Table[[Minus]Sum[
      christoffel[[i, j, k]] coords[[j]]' coords[[k]]', {j, 1, n}, {k,
        1, n}], {i, 1, n}]]
listgeodesic := 
 Table[{"d/d[Tau]" ToString[coords[[i]]'], " =", geodesic[[i]]}, {i, 
   1, n}]
TableForm[listgeodesic, TableSpacing -> {2}]

max[Tau] = 750; ivs = {0, 0, 0.088}; ics = {0, 6.5, [Pi]/2, 
  0}; r0 = 1;
computeSoln[max[Tau]i_, ivsi_, icsi_] := 
 Block[{ivs, ics, i, [Chi], tmp, soln}, ics = icsi;
  ivs = Join[{[Chi]}, ivsi];
  op1 = Table[coords[[i]] -> ics[[i]], {i, 0, n}];
  tm = metric /. op1;
  tmp = ivs.(tm.ivs); [Chi]slv = Solve[tmp == uinvar, [Chi]];
  ivs[[1]] = Last[[Chi] /. [Chi]slv];
  op = {Derivative[1][t] -> Derivative[1][t][[Tau]], 
    Derivative[1][r] -> Derivative[1][r][[Tau]], 
    Derivative[1][[Theta]] -> Derivative[1][[Theta]][[Tau]], 
    Derivative[1][[CurlyPhi]] -> Derivative[1][[CurlyPhi]][[Tau]], 
    t -> t[[Tau]], 
    r -> r[[Tau]], [Theta] -> [Theta][[Tau]], [CurlyPhi] -> 
[CurlyPhi][[Tau]]};
  deq = Table[
    coords[[i]]''[[Tau]] == Simplify[geodesic[[i]] /. op], {i, 1, 
     n}]; deq = 
   Join[deq, Table[coords[[i]]'[0] == ivs[[i]], {i, 1, n}], 
    Table[coords[[i]][0] == ics[[i]], {i, 1, n}]];
  soln = NDSolve[deq, coords, {[Tau], 0, max[Tau]i}]; soln]
uinvar = [Minus]1;
sphslnToCartsln[soln_] := 
 Block[{xs, ys, zs}, 
  xs = r[[Tau]] Sin[[Theta][[Tau]]] Cos[[CurlyPhi][[Tau]]] /. 
    soln; ys = 
   r[[Tau]] Sin[[Theta][[Tau]]] Sin[[CurlyPhi][[Tau]]] /. soln;
  zs = r[[Tau]] Cos[[Theta][[Tau]]] /. soln; {xs, ys, zs}]
udotu[solni_, [Tau]val_] := 
 Block[{x[Alpha], u[Alpha]}, 
  x[Alpha] = 
   Table[coords[[i]][[Tau]] /. solni, {i, 1, n}] // Flatten;
  u[Alpha] = D[x[Alpha], [Tau]];
  x[Alpha] = x[Alpha] /. [Tau] -> [Tau]val;
  u[Alpha] = u[Alpha] /. [Tau] -> [Tau]val;
  u[Alpha].((metric /. 
       Table[coords[[i]] -> x[Alpha][[i]], {i, 1, n}]).u[Alpha])]
coordlist[[Tau]in_] := 
 Table[ToString[coords[[i]]] <> 
   " = " <> {ToString[coords[[i]][[Tau]in] /. soln // First]}, {i, 1,
    n}]

soln = computeSoln[max[Tau], ivs, ics];

xyzsoln = sphslnToCartsln[soln];
Join[{"Final Coordinates:"}, coordlist[max[Tau]]] // TableForm
Join[{{"", "", "", "u.u values"}}, 
  Table[{"[Tau]=", ToString[i], "->", udotu[soln, i]}, {i, 0, 
    max[Tau], max[Tau]/5}]] // TableForm
{Plot[Evaluate[
   Table[coords[[i]][[Tau]] /. soln, {i, 1, n}]], {[Tau], 0, 
   max[Tau]}, AxesLabel -> {"[Tau]", "Coordinate"}, 
  PlotLegends -> {"t", "r", "[Theta]", "[CurlyPhi]"}, 
  PlotRange -> {0, 30}], 
 Show[ParametricPlot[
   Evaluate[{xyzsoln[[1]], xyzsoln[[2]]} // Flatten], {[Tau], 0, 
    max[Tau]}, AspectRatio -> 1, PlotStyle -> Gray], 
  Graphics[{Red, Circle[{0, 0}, 2]}]]}

Shape FUNCTION;

(1)/(((r)/((Sqrt[r*r0]))) - 1)

Integrate[(1)/(((r)/((Sqrt[r*r0]))) - 1), r]

max[Tau] = 2000;
ivs = {[Minus]0.08, .035, .0359};
ics = {0, 2, Pi/4, 0.2};
xyzsoln = sphslnToCartsln[computeSoln[max[Tau], ivs, ics]];
angle = ParametricPlot3D[
   Evaluate[Re[xyzsoln] // Flatten], {[Tau], 0, max[Tau]}, 
   AxesLabel -> {x, y, z}, DisplayFunction -> Identity];
sphhoriz = 
  SphericalPlot3D[{1, 2, 3, 4}, {[Theta], 0, Pi}, {[Phi], 0, 2*Pi}, 
   PlotRange -> {{-10, 10}, {-10, 10}, {-30, 30}}, 
   BoxRatios -> {1, 1, 1}, PlotTheme -> "Classic", BoxRatios -> 2, 
   Mesh -> None, PlotStyle -> Opacity[0.2]];
sphhoriz2 = 
  SphericalPlot3D[{1, 2, 3, 4}, {[Theta], 0, Pi}, {[Phi], 
    0, (3/2)*Pi}, PlotRange -> {{-10, 10}, {-10, 10}, {-30, 30}}, 
   BoxRatios -> {1, 1, 1}, PlotTheme -> "Classic", BoxRatios -> 2, 
   Mesh -> None, PlotStyle -> Opacity[0.3]];
Show[angle, sphhoriz, DisplayFunction -> $DisplayFunction, 
 PlotRange -> {{-10, 10}, {-10, 10}, {-10, 10}}]
Show[angle, sphhoriz2 , DisplayFunction -> $DisplayFunction, 
 PlotRange -> {{-10, 10}, {-10, 10}, {-10, 10}}]