Here the code for BVP, but there is a singularity how to avoid it

a:=-0.1:alpha:=0.1:
dsys:={diff(u(r),r,r)+(1-alpha)/r+r^(alpha-2)*diff(u(r),r)=a};

res1:=dsolve(dsys union {u(-1)=0,u(1)=0},numeric,output=listprocedure);

plots[odeplot](res1,[[r,(u(r))]],-5..5,axes=boxed,tickmarks = [3, 2], color = black, thickness = 1, linestyle = solid, titlefont = [Helvetica, roman, 18], labeldirections = [horizontal, vertical], labelfont = [Helvetica, roman, 24]);

here is my try

for ploting points of data ( d(n(, sum(n))

new.mw

Here is my try to integrate the expression L with trapozoid or simpson 

numerical_int.mw

Here is my code:

t := 6*Pi; lambda := .1; N := 5; M := sqrt(N(N+1))*exp(I*phi); omegap := 1

eq1:=-2*(n(phi,delta)-N)+(u(phi,delta)-M)*exp(-2*I*omegap*t/lambda)+((z(phi,delta)-conjugate(M))*exp(2*I*omegap*t/lambda)):
eq2:=-2*(1-I*delta)*u(phi,delta)+2*(n(phi,delta)-N)*exp(2*I*omegap*t/lambda)+2*M:
eq3:=-2*(1+I*delta)*z(phi,delta)+2*(n(phi,delta)-N)*exp(-2*I*omegap*t/lambda)+2*conjugate(M):
convert(eq1, rational); convert(eq2, rational); convert(eq3, rational);

f1:=evalc(eq1);f1:=evalc(eq2):f3:=evalc(eq3)
f1 := solve{f1, f2, f3};
with(plots); Impliciplot3([n(phi,delta)], delta = -10 .. 10, phi = 0 .. 2*Pi);

lambda:=0.1:N:=5:M:=sqrt(N(N+1)):omegap:=10:phi:=0:

var:={n(t),u(t)}:
dsys:={diff(n(t),t)=-2*(n(t)-N)+(u(t)-abs(M)*exp(I*phi))*exp(-2*I*omegap*t/lambda)+conjugate((u(t)-abs(M)*exp(I*phi))*exp(-2*I*omegap*t/lambda)),diff(u(t),t)=-2*(1-I*delta)*u(t)+2*(n(t)-N)*exp(-2*I*omegap*t/lambda)+2*abs(M)*exp(I*phi)}:
res1:=dsolve(dsys union {n(0)=0,u(0)=0},numeric,output=listprocedure):

# (this command need to be changed) P1:=plots[odeplot3](res1,[[t,(n(t))]],0..10,axes=boxed,tickmarks = [3, 2], color = black, thickness = 4, linestyle = solid, titlefont = [Helvetica, roman, 18], labeldirections = [horizontal, vertical], labelfont = [Helvetica, roman, 24]);