@toandhsp See the Reply by Markiyan. The equation has two real solutions. Using RealDomain, the equation has one real solution. In my opinion, RealDomain is correct to reject the solution -1-sqrt(2). If RealDomain were simply a filter to discard non-real solutions, it would be rather trivial.

@dcanright Matrices are not just intended for algebra. This is just a little bug that prevents you from using a one-dimensional initializer. Using a two-dimensional initializer of strings works fine:

Matrix(2,2, [["A", "B"], ["C", "D"]]);

@Markiyan Hirnyk In the 3d plot that you show, one function is linear and the other is piecewise linear. To me, those are not believable representations of the functions.

@radzys Unfortunately, using your suggested fifth boundary condition (D@@3)(phi)(L) = (D@@2)(phi)(L) leads to the error message "Matrix is singular", probably because it leaves the possibility of an identically zero solution.

@Markiyan Hirnyk 

The collect(combine(expand(...))) does nothing is this case. Indeed, they leave p unchanged. Only the coeff(p, x^n)  (rather than the OP's coeff(p, x, n)) is needed. The following also work:

coeff(p, x^(n+2));

coeff(coeff(p, x^(n+2)), y^(n+2));

Your system contains parameters omega and k_r for which no numerical values are given. Do you have values for these?

Runge-Kutta or rkf45 can only be used for initial value problems (IVPs). This is a boundary value problem (BVP). Better to just let dsolve choose the method.

You have numerous syntax errors which I'll show you how to correct once you supply the values of the parameters.

@wo0olf 

Use the output of dsolve as a function:

sol1(.0001);

@hacker9130368 Going from intensity=9000000 to intensity=10000000, I do get an increase in b and a decrease in c.

With intensity=9000000:

With intensity=10000000:

If you got no difference, then you probably need to set Digits to a higher value. I used Digits:= 15.

How did you embed the animation in your post? Is it a .GIF file? It is unusual in that the animation only goes through one cycle.

@mattcanderson1 Regarding Maple's use of the GPU, see ?CUDA in the Maple help system. However, this is unlikely to be any help in this thread's situation.

@manueleleonelli 

How about this?:

lower:= x^2*y, x= 0..1, y= 0..1, grid= [10$2], filled:
upper:= x^2*y-1, x= 0..1, y= 0..1, grid= [10$2], filled:
plots:-display(
     [plot3d(lower, style= wireframe),
      plot3d(lower, style= patchnogrid),
      plottools:-translate(              
           plot3d(upper, style= patchnogrid, color= pink),
           0, 0, 1
      )],
     orientation= [130,65], axes= none,
     view= [(0..1) $ 3], transparency= 0.2
);

@acer Using a= 1/2, b= 1/2, A= 0.1, B= 0.9, I get isconvex returning false without resorting to the definition. Note that f has singularities at x=0 and x=1, so using A=0, B=1 is questionable.

The false result is wrong: the function is convex.

Using eliminate instead of solve and ignoring the parametric aspect, we get the following weird result:

eliminate({
     v*(1-a) = alpha[3], a*u*v*b = alpha[2], v*a*u*(1-b) = alpha[1],
     (1/2*(-1+a))*a*p*v = alpha[33], (-1+a)*(-1+b)*a*p*v*u = alpha[13],
     (1-a)*a*p*v*u*b = alpha[23], (1/2)*a*u*v*b*(-b*u*p+b*a*u*p+b*g-g) = alpha[22],
     -a*u*v*b*(p*u-u*p*a-g-b*u*p+b*a*u*p+b*g) = alpha[12],
     (1/2)*a*u*v*(-1+b)*(-b*u*p+b*a*u*p+b*g+p*u-u*p*a) = alpha[11]
     }, {a, b, g, p, u, v}
);

@hacker9130368 Thanks for the great help on converting help files! Ironically, there is no normal help page for HelpTools:-Database:-ConvertAll. I got it to work though. After using it, you need to close your help window and re-open it; I guess it needs that to load the new index.

@mkj256 I it's more likely that it doesn't like fractions, which is how the first column is represented, To check, multiply by 1.0.