Where did you read or hear about a package named Lamp? It could be a package made by a user of Maple. It is not a package that is distributed with Maple.

You cannot use brackets [ ] instead of parentheses in Maple. You have several brackets that should either simply be removed or, if needed, replaced by parentheses.
An example:
You have

f25 :=r->[diff(lambda(r), r)*r-lambda(r)-(n^2+2+alpha^2*r^2)*mu(r)+omega^2*rho*r^2];

You should write this simply without the brackets:

f25:=r->diff(lambda(r), r)*r-lambda(r)-(n^2+2+alpha^2*r^2)*mu(r)+omega^2*rho*r^2;

Another thing: Are you sure that you need to define f25 (as an example) as a function (procedure in Maple)?
Wouldn't this do the job:

f25:=diff(lambda(r), r)*r-lambda(r)-(n^2+2+alpha^2*r^2)*mu(r)+omega^2*rho*r^2;

Then when used write f25 instead of f25(r).
The problem with the definition of f25 as a procedure is that f25(1.2345) wouldn't make any sense since that call to f25 would try this
diff(lambda(1.2345), 1.2345)
which obviously is nonsense.

@matlar If you found a solution with those initial conditions it couldn't have been on all of t = 0..100.
The solutions I found above are all numerical. Both dsolve/numeric and DEtools[DEplot] use per default the method rkf45.
If you want to use RK4 you can give dsolve the optional arguments method=classical[rk4] and stepsize=0.01 (or whatever stepsize you like).

@carriewong You have 'local' twice. Remove the first occurrence of 'local'.
You should include evals and evecs among the locals by replacing the colon after Dg with a comma.

 

When I enter this command in Maple
exp(2*t);
 I get the following 

In other words, I don't see any parentheses.
I see that you are using Maple 13. I also tried in Maple 12 where the output looked quite the same.

@oceanxinlie The transformation of your boundary value problem:
ode:=diff(y(x),x,x)=a*abs(y(x))^q*signum(y(x))+b*x;
with y(0) = 0, D(y)(L) = 0 can be transformed to the form:

diff(v(z),z,z)=abs(v(z))^q*signum(v(z))+B*z;
with v(0) = 0, D(v)(1) = 0.

and so can the ode where signum(y(x)) is replaced by -1, but the version considered here is the one that corresponds to your piecewise expression.
The code:

restart;
ode:=diff(y(x),x,x)=a*abs(y(x))^q*signum(y(x))+b*x;
PDEtools:-dchange({y(x)=A*v(z),x=L*z},ode,[v,z]);
ode1:=expand(op(solve(%,{diff(v(z),z,z)}))) assuming A>0;
ode1:=simplify(eval(ode1,A=(L^2*a)^(1/(1-q)))) assuming L>0,a>0,b>0,q>0,q<1;
# Your parameters are:
a:=0.00003019: b:=9.542*10^(-13): L:=2945: q:=0.337:
## So the two odes are
ode;
ode1;

@oceanxinlie Your last  question:

"you have mentioned the ode:

ode:=diff(y(x),x,x)=-abs(y(x))^(1/3)+x;

Can you solve it by code in maple?  Is it using the simple shooting method?  I have a lot of places confused. "

Yes, a simple shooting method is used and I already gave the full code in Maple in the reply with the title "Not quite yet" at the end. Notice that the code starts with a restart command. That is a clear indication that everything is in there, nothing is left out except a detailed explanation of the code. I didn't care to do that since I hadn't (and still haven't) been successful in using it on your system transformed to a similar form. That transformation I shall give in a second reply.

@oceanxinlie The ode given by

ode:=diff(y(x),x,x)=abs(y(x))^0.337*signum(y(x))+b*x;

  is the same as the one you mention (after the transformation I mentioned above)

f''(x)=piecewise(f(x)>=0, 3.019*(10^(-5))*f(x)^0.337, [-3.019*(10^(-5))*(-f(x))^0.337])+9.541*(10^(-13))*x

because
abs(y(x))^0.337*signum(y(x)) = y(x)^0.337 for y(x) >= 0
and
abs(y(x))^0.337*signum(y(x)) = -(-y(x))^0.337 for y(x) < 0.

@oceanxinlie Solving the problem

ode:=diff(y(x),x,x)=abs(y(x))^0.337*signum(y(x))+b*x;

considered on x = 0..1 using the shooting method shown above begins to run into problems already at values of b less than 0.1.
The b you have (once your ode is transformed as mentioned above) is 0.549269696597376e-5. For that value of b the graph of  y(0) as a function of y(1) shows what looks like a discontinuity where the sign changes:

It is probably not really a discontinuity, but "only" a numerical problem. It is not surprising, however, that the derivative of y(0) w.r.t. to the parameter y(0)=y0 seems to be infinite at a point where y = 0 because the exponent of y(x) in the ode is 0.337 < 1.
Surely it doesn't help that at x = 1 the value of y(1) = y0 must be chosen close to zero at about -2.5*10^(-16).
Shown below are the graphs of two "solutions" found by shooting. The value of y(1) = y0 differs by only 1e-30 between the two:

It is remarkable how bad it is.
The problem here described can be seen for much higher values of b as mentioned above.

@oceanxinlie What made me think that I was on the right track was the success for this simplified version, where the ode is simply

ode:=diff(y(x),x,x)=-abs(y(x))^(1/3)+x;

and will be considered on the interval x = 0..1.
Since we only accept y(x)<=0 this abs version should be fine.
Your problem is

ode:=diff(y(x),x,x)=-a*abs(y(x))^0.337+b*x;

on x = 0..2945 with a and b as given by you and boundary conditions y(0)=0, D(y)(2945) = 0.
Your problem can by a change of variables y(x) = A*v(z), x=L*z, be transformed into
 

ode1:=diff(v(z),z,z)=-abs(v(z))^0.337+B*z;

on the interval z = 0..1 and with v(0)=0, D(v)(1)=0, where B = 0.549269696597376e-5.
Thus your problem is close to the my toy problem except for the coefficient B (and the unimportant difference between 0.337 and 1/3). That B is so small seems to be a problem. Right now I don't know why, but will continue looking at it for a while.
##
Here follows the code for the toy problem. It uses a simple shooting method and shoots from x = 1 trying to make y(0) = 0 by varying the value of y(1) = y1. It uses fsolve with a guessed value for y1 found by looking at the animation provided by the code.
dsolve with the parameters option is used.
 

restart;
ode:=diff(y(x),x,x)=-abs(y(x))^(1/3)+x;
plots:-shadebetween(-x^3,0,x=0..1,color=gray); bg:=%:
res:=dsolve({ode,y(1)=y0,D(y)(1)=0},numeric,parameters=[y0],output=listprocedure);
Y:=subs(res,y(x));
Q:=proc(y0,b::truefalse:=true) if not y0::realcons then return 'procname(_passed)' end if;
  res(parameters=[y0]); 
  if b then Y(0) else res end if
end proc;  
plots:-odeplot(Q(-0.1,false),0..1); #Testing Q
plots:-animate(plots:-odeplot,['Q(y0,false)',[x,y(x)],0..1],y0=-0.2..-0.1,background=bg);
fsolve(Q,-0.116); #Finds the value of y(1) for which y(0)=0.
res(parameters=[%]); #Setting the parameter y1.
plots:-odeplot(res,[x,y(x)],0..1,caption="The solution");

@oceanxinlie I'm pretty sure I now have a solution. I'll be gone for a few hours, but will report back then.

Use VectorCalculus:-Jacobian as in
 

VectorCalculus:-Jacobian(L1, L2);
J:=unapply(%,L2);

where L1 is the list of right hand sides and L2 is the list of names (variables).
Follow up with the unapply line and then J will be a function of the variables.

You need pdsolve instead of solve.

@oceanxinlie I didn't use the expression for u that you use above.
I used

u:=piecewise(y(x)>=0,y(x)^0.337,0);

and remarked at the bottom the note: 
"This 'desperate' attempt to get a positive solution by rewriting the ode ends up just giving the same as the one where u = 0."
So that means I discarded that idea!
Now your expression u is
 

u:=piecewise(y(x)<0,0,-(-y(x))^0.337);

but that won't work because (-y(x))^0.337  will give imaginary results for y(x) > 0.
Your previous expression

piecewise(y(x)>=0,y(x)^0.337,-((-y(x))^0.337))

made more sense.
When I asked if y(x) can change sign, I didn't mean to pose a mathematical question, but rater this: Physically (apart from the actual mathematical mode) can y(x) change sign?
It appears that you actually answer 'No' when you say:
"Now I knew y(x)<0,  so the new ode is shown below:

y''(x)+0.00003019*[-y(x))]^0.337=[9.541*10^(-13)]*x  "

But (-y(x))^0.337 = abs(y(x))^0.337 for y(x) <=0, so in principle you might as well use
 

ode1:=diff(y(x),x,x)+0.00003019*abs(y(x))^0.337=9.542*10^(-13)*x; 

which, however, will give you a positive solution y(x), which is actually roughly the same as -y(x) from

ode2:=diff(y(x),x,x)-0.00003019*abs(y(x))^0.337=9.542*10^(-13)*x; 

except for the sign. ode2 gives a negative solution as seen in my answer.

If y(x) solves ode1 then z(x) = -y(x) will solve
-diff(z(x),x,x)+0.00003019*abs(z(x))^0.337=9.542*10^(-13)*x;
i.e.
diff(z(x),x,x)-0.00003019*abs(z(x))^0.337=-9.542*10^(-13)*x;

The right hand sides have different signs, but they are small in absolute value.
The maximal relative difference of the solutions of ode1 and ode2 is 1.7e-5.

@oceanxinlie Is it expected that y(x) actually changes sign?
The immediate problem encountered when trying dsolve on your revised ode is that the piecewise expression isn't differentiable at zero:

ode:=diff(y(x),x,x)-0.00003019*piecewise(y(x)>=0,y(x)^0.337,-((-y(x))^0.337))=9.541*10^(-13)*x;
dsolve({ode,bcs},numeric,method=bvp[midrich]); # ERROR

plot(piecewise(y>=0,y^0.337,-(-y)^0.337),y=-1..1);
diff(piecewise(y>=0,y^0.337,-(-y)^0.337),y);