Some part of the Circular path is not visible when view sidewise. But when view from top the whole circular path is visible in my 3D plot matplotlib

Question

Using matplotlib I plotted the circular path around a sphere. Some parts of the plot on the right and left sides are visible when viewed from the sides. But when viewed from the top the whole circular path is visible as in image 2. Someone, please check what is wrong with the code.

The code used is as under:

from scipy.integrate import solve_ivp
import numpy as np
from math import sin, cos,pi,sqrt

μ = 379312077e8           # m^3/s^2 # G= 6.67408e-11m^3/kg*s^ # M = 5.68e26 # μ=GM
R = 60268e3               # metre
g_10 = 21141e-9
Ω = 17e-5                 # rad/s = 9.74e-7   

# charge by mass ratio
ρ = 1e3                   # 1gm/cm^3 Headmen 2021= 10^3 kg/m^3, 
b = 1e-11                  # 1nm = 1e-3 micro
V = 10                    # Volt
ε = 8.85e-12              # Farad/metre        
m = (4/3)*pi*(b**3)*ρ     #
q = b*V*4*pi*ε            # C
β1 = q/m
β  = (3*ε*V)/(ρ*b**2)
L = β*((g_10*R**3*Ω)/μ)

def odes(t,p):
# assigning each ODE to a vector element
    r,x,θ,y,ϕ,z = p

# constants
    μ = 379312077e8                 # m^3/s^2,G=6.67408e-11m^3/kg*s^ # M = 5.68e26 #kg 
    R = 60268e3                     # metre
    j_2 = 1.629071e-2
    g_10 = 21141e-9                      #
    Ω = 1.7e-4                           # rad/second
    μ_0 = 4*pi*1e-7 
    B_θ = μ_0*(R/r)**3*g_10*sin(θ)
    B_r = μ_0*2*(R/r)**3*g_10*cos(θ)
    β = (3*ε*V)/(ρ*b**2)                # q/m = -3.46 x 103 C/kg, 

# defining the ODEs
    drdt = x
    dxdt = r*((y**2 +((z+Ω)**2)*(sin(θ)**2))-β*z*sin(θ)*B_θ)-(μ/r**2)*(1-(3/2)*j_2*((R/r)**2)*((3*cos(θ)**2)-1))
    dθdt = y
    dydt = (-2*x*y + r*((z+Ω)**2)*sin(θ)*cos(θ)+ (β*r*z*sin(θ)*B_r))/r + (3*μ/r**2)*j_2*((R/r)**2)*sin(θ)*cos(θ)
    dϕdt = z
    dzdt = (-2*(z+Ω)*(x*sin(θ)+r*y*cos(θ)) + (β*(x*B_θ-r*y*B_r)))/(r*sin(θ))

    return np.array([drdt,dxdt,dθdt,dydt,dϕdt,dzdt])

# time window
t_span = (0,86400*30)
t = np.linspace(t_span[0], t_span[1], 1000) 
    
#initial conditions
r0 = 1.12*R
y_0 = (μ/r0)**0.5
y_w = (y_0/(2*pi*r0))*360
p0 = np.array([1.12*R,0.0,90.0*(pi/180),0.0*(pi/180), 0.0*(pi/180), y_w*(pi/180)])

sol = solve_ivp(odes,t_span,p0,method="RK45",t_eval=np.linspace(t_span[0],t_span[1], 1000))

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

# 3D Cartesian
fig = plt.figure()
x1= sol.y[0,:] * np.sin(sol.y[2,:])* np.cos(sol.y[4,:])
y1= sol.y[0,:] * np.sin(sol.y[2,:])* np.sin(sol.y[4,:])
z1= sol.y[0,:] * np.cos(sol.y[2,:])
ax = plt.axes(projection="3d")
ax.plot3D (x1,y1,z1, color="red", linewidth ='0.5')

np.sin(θ)* np.cos(ϕ) 
ax.set_xlabel('X-Axis')     # plt.xlim(-1e3,7e7)     
ax.set_ylabel('Y-axis')#plt.ylim(-1e4,7e8)
ax.set_zlabel('Z-axis')     # ax.set_zlim(-1e3,1e5)

phi = np.linspace(0, 2*np.pi, 100)     
theta = np.linspace(0, np.pi, 100)
xm = R*np.outer(np.cos(phi), np.sin(theta))
ym = R*np.outer(np.sin(phi), np.sin(theta))
zm = R*np.outer(np.ones(np.size(phi)), np.cos(theta))
plot = ax.plot_surface( xm, ym, zm, rstride=1, cstride=1, cmap=plt.get_cmap('jet'), 
                   linewidth = 0.5, antialiased = False, alpha = 1)
plt.show()