kill(all)$globalsolve:true$

assume(r>0,R>0,L>0);

a:ω**2*x;
g:K/x**2;

X:ρ*(a-g);

N:integrate(X,x,r,R+L);

diff(N,r);

solve(diff(N,r)=0,r);

N_0:ev(N,r:R);

N_0:ratsimp(N_0);

sol:solve(N_0=0,L);

L_0:rhs(sol[2]);

ev(L_0/R,R:6378e3,K:9.81*R**2,ω:2*π/24/3600,eval,numer);

plot2d(ev(N,ρ:1550,R:6378e3,L:1.442905955394484e8,K:9.81*R**2,ω:2*π/24/3600,eval,numer),[r,6378000,6378000*23.62])$