@_Maxim_ Just a comment about the last point. Notice the difference between these two calls, where the sequence defined as a is evaluated by evalf at the setting of Digits, and the second where evalf is given the explicit sequence Pi,7 whereby Pi is evaluated at Digits=7:

restart;
Digits:=10:
a:=Pi,7;
evalf(a);
evalf(Pi,7);

@tomleslie There is no doubt that the system dsys3 satisfies the requirements for existence and uniqueness.
Solving for the second derivatives gives right hand sides that are just a sum of terms including besides sin(t) and constants only polynomials of order 3 in the dependent variables and their derivatives. 

But how large the maximal interval of definition of the solution is, is certainly not clear. It may be (1) finite, and (2) very small. Thus it might not even include  t = 1.

I notice that your equations eq3 and eq4 contains integrals. This one is from eq4:
int(cosh(4.6941*x)-cos(4.6941*x)-.9825251551*sinh(4.6941*x)+.9825251551*sin(4.6941*x), x = 0 .. .999)

and it just evaluates to .8369623967. The similar integral found in eq3 just evaluates to .2866222667.

Is this as intended?

@Carl Love Since the values in the output = array(...) are the values of the independent variable, it doesn't make much sense to ask for t=1 three times.

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.