e=0.1; (* mimośród *)
(Subscript[ϕ, 0])=0; (* początkowa szerokość geograficzna *)
(Subscript[λ, 0])=0; (* początkowa długość geograficzna *)
α=1; (* początkowy kierunek *)
(Subscript[t, max])=100; (* czas maksymalny *)

(Subscript[θ, 0])=ArcTan[(1-(e^2))Tan[(Subscript[ϕ, 0])]];
s=NDSolve[{(((Tan[ϕ[t]])^2)+1-(e^2))ϕ''[t]+(e^2)Tan[ϕ[t]]((ϕ'[t])^2)+Tan[ϕ[t]]((λ'[t])^2)==0,λ''[t]==2Tan[ϕ[t]]ϕ'[t]λ'[t],ϕ[0]==(Subscript[θ, 0]),λ[0]==(Subscript[λ, 0]),ϕ'[0]==((Cos[α])/(Sqrt[((Sin[(Subscript[θ, 0])])^2)+((Cos[(Subscript[θ, 0])])^2)(1-(e^2))])),λ'[0]==((Sin[α])/(Cos[(Subscript[θ, 0])]))},{ϕ,λ},{t,0,(Subscript[t, max])}];Show[ParametricPlot3D[Table[{Cos[ϕ]Cos[λ],Cos[ϕ]Sin[λ],(Sqrt[1-(e^2)])Sin[ϕ]},{ϕ,-(π/2),(π/2),(π/6)}],{λ,0,2π},PlotStyle->Gray],ParametricPlot3D[Table[{Cos[ϕ]Cos[λ],Cos[ϕ]Sin[λ],(Sqrt[1-(e^2)])Sin[ϕ]},{λ,-π,π,(π/6)}],{ϕ,-(π/2),(π/2)},PlotStyle->Gray],ParametricPlot3D[{Cos[ϕ[t]]Cos[λ[t]],Cos[ϕ[t]]Sin[λ[t]],(Sqrt[1-(e^2)])Sin[ϕ[t]]}/.s,{t,0,(Subscript[t, max])},PlotStyle->Black],Axes->False,Boxed->False]