@Carl Love 

Clear explanation! Thank you so much, I think I have to come back to do more work on the validation of this ODEs and BCs. 

But I did learn a lot here, thank you again!

@Carl Love 

I am sorry, what do you mean by force the value of the constants?

You mean the extra BC will induce conflict with these constants?

If I omit one of the BCs, then

The answer will be like:

But as I know this system should be somehow symmetrical to the centre ( x=50 ), I don't know what's wrong here...

Any suggestion?

@Preben Alsholm 

Thank you! 

But why I still can't plot it with the actual boundary conditions, is there anything wrong?

Following is the codes：

restart; 
sys_ode := 2*C__5*(diff(w(x), x))-C__1*(diff(u(x), x$2))-2*C__4*(diff(w(x), x$3)) = (1/2)*q*(x-50)^2-1250*q, 2*C__2*u(x)-C__1*(diff(w(x), x$3))-2*C__3*(diff(u(x), x$2)) = 0;
ics := w(0) = 0, w(100) = 0, (D(u))(0) = 0, ((D@@2)(w))(0) = 0, (D(u))(100) = 0, ((D@@2)(w))(100) = 0;
map(int, sys_ode[1], x = 0 .. x, continuous);
ode1 := eval(%, {ics});
ode2 := sys_ode[2]；
res := dsolve({ics, ode1, ode2})：
U, W := op(subs(res, [u(x), w(x)]))：
Digits := 200；
plot(eval(W, {C__1 = 0.12e9, C__2 = 0.2e10, C__3 = 0.2e9, C__4 = 0.3e8, C__5 = 0.3e8, q = 1}), x = 0 .. 100, size = [2000, default])

@Carl Love 

@Preben Alsholm 

In the original question, I used:

ics := (D(u))(0) = 1, ((D@@2)(w))(0) = 0, (D(u))(100) = 1, ((D@@2)(w))(100) = 2                     (FOUR there ) 

But I just used these values 1,0,1,2 to try if I can get the solution, the actural boundary conditions are:

ics :=w(0) = 0, w(100) = 0, (D(u))(0) = 0, ((D@@2)(w))(0) = 0, (D(u))(100) = 0, ((D@@2)(w))(100) = 0                         (SIX there)

These are got from the physical meaning, which:

 ---w(0) = 0, w(0) = 0: 

Simply supportted beam, hence, no deflection in the supports

---(D(u))(0) = 0, ((D@@2)(w))(0) = 0, (D(u))(100) = 0, ((D@@2)(w))(100) = 0

From:

and 

Note:

 and  can be arbitrary, soand  is 0.

Further: M(0)=0, and M(l)=0, where the l could be 100, which is used in your solutions.

So the question is,

I used the actual boundary condition, but it seems the outcomes are not reasonable, is it possible I used the wrong boundary condition?

@Preben Alsholm 

I changed the Function Q (x) to  (1/2)*q*(x-50)^2-1250*q, and adjusted the boundary conditions to a more reasonable set,  it seems not working now ... Can you tell me why?

restart;
sys_ode := 2*C__5*(diff(w(x), x))-C__1*(diff(u(x), x$2))-2*C__4*(diff(w(x), x$3)) = (1/2)*q*(x-50)^2-1250*q, 2*C__2*u(x)-C__1*(diff(w(x), x$3))-2*C__3*(diff(u(x), x$2)) = 0;
ics := w(0) = 0, (D(u))(0) = 0, ((D@@2)(w))(0) = 0, (D(u))(100) = 0, ((D@@2)(w))(100) = 0; 
map(int, sys_ode[1], x = 0 .. x, continuous); 
ode1 := eval(%, {ics}); 
ode2 := sys_ode[2]; 
res := dsolve({ics, ode1, ode2}); 
U, W := op(subs(res, [u(x), w(x)])); 
Digits := 200; 
plot(eval(U, {C__1 = 0.12e9, C__2 = 0.2e10, C__3 = 0.2e9, C__4 = 0, C__5 = 0.3e8, Q = 1}), x = 0 .. 100, size = [2000, default])

@Carl Love 

Truly thanks for your help!

But I changed the boundary condition to a more complex one, the codes seem not working?

Changed_Boundary_Explore_BVP.mw

By the way,  is it applicable if the Q(x) is a piecewise function?

@Carl Love 

It really helps!

@Carl Love 

Yes, you are right, I missed this one, w(0) = 0, which is indeed the boundary my model should fulfil. 

Also, yes, I forgot that I have limited the x in 0 .. 100, but the length I actually need is around 5000, at this order of magnitude roughly. 

Sorry for the late reply.

@Carl Love 

The interval of variable ' x ' is from ' 0 ' to ' 5000 '. 

Thank you again! The complex solutions might not be expected, they will not that easy to explain with physical meanings...

@Carl Love 

 I love Maple!

Thank you again!

Following is the acceptable ranges:

C__1 :  1.2 ~ 1.8 *(10^8)

C__2 :  2 ~ 10 *(10^9)

C__3 :  2 ~ 6 *(10^8)

C__4 :  0 ~ 5 *(10^7)

C__5 :  3 ~ 10 *(10^7)

@Carl Love 

Thanks for your helpful response, with sharp insight!

Yes, I was trying to make Maple give me the answer with any arbitrary Q(x) ... 

Is it unrealistic?? haha ... sorry I am new to Maple...

And I do really want to know how these C_s (coefficients, which is related to my physical model) affect the final outcomes.

Do you have any suggestion for the numeric solver? or any sample? just for this question, it's too time-consuming to learn all the operations...