systeme={r'[t]==dr[t],
dr'[t]==-Sin[theta[t]]^2/r[t]^2*dphi[t]+r[t]*dphi[t]^2*Sin[theta[t]]^2+r[t]*dtheta[t]^2,
theta'[t]==dtheta[t],
dtheta'[t]==(2*Cos[theta[t]]*Sin[theta[t]]*dphi[t]/r[t]^2+r[t]*dphi[t]^2*Sin[theta[t]]
            *Cos[theta[t]]-2*dr[t]*dtheta[t])/r[t],
phi'[t]==dphi[t],
dphi'[t]==((dr[t]*Sin[theta[t]]-r[t]*dtheta[t]*2*Cos[theta[t]])/r[t]^3-2*r[t]*dphi[t]
            *dtheta[t]*Cos[theta[t]]-2*dr[t]*dphi[t]*Sin[theta[t]])/(r[t]*Sin[theta[t]])}
            
initial={r[0]==1,dr[0]==0,theta[0]==Pi/2,dtheta[0]==0,phi[0]==0,dphi[0]==0.1}
tmax=200
solution=NDSolve[Join[systeme,initial],{r,dr,theta,dtheta,phi,dphi},{t,0,tmax},
                Method->{"ExplicitRungeKutta",DifferenceOrder->5},
                MaxSteps->Infinity,AccuracyGoal->6]
            
Plot[r[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,2}},AxesLabel->{"t","r"}]
Plot[theta[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,Pi}},AxesLabel->{"t","theta"}]
Plot[phi[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,2*Pi}},AxesLabel->{"t","phi"}]
Plot[Sqrt[dr[t]^2+(r[t]*dtheta[t])^2+(r[t]*Sin[theta[t]]*dphi[t])^2]/.solution,{t,0,tmax},
                PlotRange->{{0,tmax},{0,0.2}},AxesLabel->{"t","V"}]
ParametricPlot3D[{r[t]*Sin[theta[t]]*Cos[phi[t]],r[t]*Sin[theta[t]]*Sin[phi[t]],
                    r[t]*Cos[theta[t]]}/.solution,{t,0,tmax},PlotRange->{{-2,2},{-2,2},{-2,2}},
                    AxesLabel->{"x","y","z"},PlotPoints->1000]
initial={r[0]==1,dr[0]==0,theta[0]==Pi/2,dtheta[0]==0.02,phi[0]==0,dphi[0]==0.02}
tmax=1000
solution=NDSolve[Join[systeme,initial],{r,dr,theta,dtheta,phi,dphi},{t,0,tmax},
                Method->{"ExplicitRungeKutta",DifferenceOrder->5},
                MaxSteps->Infinity,AccuracyGoal->6]
            
Plot[r[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,2}},AxesLabel->{"t","r"}]
Plot[theta[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,Pi}},AxesLabel->{"t","theta"}]
Plot[phi[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,Pi}},AxesLabel->{"t","phi"}]
Plot[Sqrt[dr[t]^2+(r[t]*dtheta[t])^2+(r[t]*Sin[theta[t]]*dphi[t])^2]/.solution,{t,0,tmax},
                    PlotRange->{{0,tmax},{0,0.1}},AxesLabel->{"t","V"}]
ParametricPlot3D[{r[t]*Sin[theta[t]]*Cos[phi[t]],r[t]*Sin[theta[t]]*Sin[phi[t]],
                    r[t]*Cos[theta[t]]}/.solution,{t,0,tmax},PlotRange->{{-2,2},{-2,2},{-2,2}},
                    AxesLabel->{"x","y","z"},PlotPoints->1000]
initial={r[0]==1,dr[0]==0,theta[0]==Pi/2,dtheta[0]==0.05,phi[0]==0,dphi[0]==0.02}
tmax=1000
solution=NDSolve[Join[systeme,initial],{r,dr,theta,dtheta,phi,dphi},{t,0,tmax},
                Method->{"ExplicitRungeKutta",DifferenceOrder->5},
                MaxSteps->Infinity,AccuracyGoal->6]
            
Plot[r[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,2}},AxesLabel->{"t","r"}]
Plot[theta[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,Pi}},AxesLabel->{"t","theta"}]
Plot[phi[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,Pi}},AxesLabel->{"t","phi"}]
Plot[Sqrt[dr[t]^2+(r[t]*dtheta[t])^2+(r[t]*Sin[theta[t]]*dphi[t])^2]/.solution,{t,0,tmax},
                    PlotRange->{{0,tmax},{0,0.1}},AxesLabel->{"t","V"}]
ParametricPlot3D[{r[t]*Sin[theta[t]]*Cos[phi[t]],r[t]*Sin[theta[t]]*Sin[phi[t]],
                        r[t]*Cos[theta[t]]}/.solution,{t,0,tmax},PlotRange->{{-2,2},{-2,2},{-2,2}},
                        AxesLabel->{"x","y","z"},PlotPoints->1000]
