********************************
Orbits

Clear["Global`*"]

power=0.7

values={m->1, el->1, k->2,energy->-1};

Veff=((1/2) el^2/(m r^2) + (-k/r^power))

Vkepler= (-k/r^power);
Vcent=(1/2) el^2/(m r^2);
Veff=Vkepler + Vcent

Plot[{Vkepler,Vcent, Veff,energy} /.values //Evaluate, {r,0.01,3}
,AxesOrigin->{0,0},
PlotStyle->{Red,Green,Blue,Black},
PlotRange->Automatic,
PlotLabel->"Energy must be somewhat above blue line!"

]

T=energy-Veff

eq1= (1/2) m r'[t]^2 - T ==0

dr2=r'[t]^2 /.Solve[eq1,r'[t]][[1]] //Apart

Plot[dr2 /.values,{r,0.1,4}];

root1=FindRoot[0==dr2 /.values //Evaluate,{r,0.02}]

root2=FindRoot[0==dr2 /.values //Evaluate,{r,1}]

{rmax,rmin}=  {root1,root2} //Map[Last,#,{2}]& //Flatten

epilog={Thickness[0.005],Red,Circle[{0,0},r] /.root1,Circle[{0,0},r] /.root2};
epilog //Graphics

r0=1;
theta0=0;
t0=0;
timeStep=0.01;
rSign=+1;
xarray={};
print=False;

Do[dr0=(rSign Sqrt[dr2] timeStep)//.values//.r->r0;
r1=r0+dr0;
dtheta0=el/(m r0^2)//.values;
theta1=theta0+dtheta0 timeStep;
test=dr2//.values//.r->r1;
(*If[test≤0,r1=r0-1 dr0;rSign=-rSign;];*)If[test≤0,r1=r0-1 dr0;rSign=-rSign;];
t0=t0+timeStep;
r0=r1;
theta0=theta1;
xarray=Append[xarray,{r1 Cos[theta1],r1 Sin[theta1]}];
If[print,Print[t0,"   ",r0,"   ",dr0,"   ",theta0];],{i,1,(10^4)}]

ListPlot[xarray,Joined->True,Epilog->epilog,PlotRange->All,AspectRatio->1]