@Aaeru Michi By default dsolve/numeric/bvp finds an approximate solution itself. That one is based entirely on the boundary conditions. In your case with

R(Pi) = 1/500^2, r(Pi) = 500, D(R)(Pi) = 0, D(r)(Pi) = 0, A(0) = 0;

it uses the approximate solution

approxsoln=[A(x) = 0, R(x) = 1/250000, r(x) = 500]

As you will notice this approximate solution just consists of constant functions, but they certainly satisfy the boundary conditions. Takíng the ode system itself into account you may be able to do better by explicitly giving an option of the form

approxsoln = [A(x)= Ac(x), R(x)=Rc(x), r(x)=rc(x)]

where each of Ac, Rc, and rc are concrete functions of x that you think might get the numerical solution started well.
Although I have tried I haven't been successful with my guesses for Ac, Rc, or rc.

I tried solving for r(x) using the ode with diff(r(x),x,x) and just setting R(x) = 1/500^2. Specifically I did this:
 

ode_r:=diff(r(x), x, x) = -(1/2)*r(x)^3*R(x)+r(x)-diff(r(x), x)*cos(x)/sin(x)+diff(r(x), x)^2/r(x);
res:=dsolve({eval(ode_r,R(x)=1/500^2),r(Pi)=500,D(r)(Pi)=0},numeric,output=listprocedure);
res(0); #Problem very, very close to 0.
rp:=subs(res,r(x));
evalf(Int(rp,0..Pi));
plot(Int(rp,0..x),x=0..Pi,numpoints=10,adaptive=false);
plots:-odeplot(res,[x,r(x)],0..Pi);

That was done to get an idea of a guess for r(x) and A(x).  I haven't had any successes so far, though.

@Markiyan Hirnyk You are quite right. But if eager0626 wants to see the roundoff then plotting the difference makes sense.

@torabi You have a system which is a variant of the famous Lorenz equations. But you ought to explain for what values of t you measure y(t) for different values of A.
In the meantime here is a nice picture to look at. A is 25.

It is generated by the following using your code except for changing A:=1 to A:=25.
 

plots:-odeplot(dsol5,[x(t),y(t),z(t)],10..110,numpoints=10000,caption="A = 25");

@Christopher2222 

You can try this as an example:
 

restart;
type(a[1],symbol);
type(a[1],name);

@John Fredsted Another simple way using that exp(1/i) > 1 for all i >=1 (thus P is an increasing function of n):

P := product(exp(1/i),i = 1..n);
fsolve(P=100,n);

The output from fsolve is 55.64520627, so the answer is 56.

@_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;