@acer I tried evalindets in Maple 2016 and also in Maple 2015.2:

restart;
A := [ [2,3], [4,5] ,[6,7]];
#indets(A,list(realcons));
evalindets(A,list(integer),(`/`)@op); #OK
evalindets(A,list(realcons),(`/`)@op); # Surprise: 35/36
#That result seems to be the result of recursing as in:
evalindets(A,list(integer),(`/`)@op);
evalindets(%,list(realcons),(`/`)@op);
###
## Doing the same in Maple 18.01:
evalindets(A,list(realcons),(`/`)@op); # [2/3,4/5,6/7]

### Were the changes to evalindets made in Maple 2015 intended to have this recursive effect or is this a bug?

With the same A and an unassigned f:
evalindets(A,list(algebraic),f);
## In Maple 2016 and 2015:  f([f([2, 3]), f([4, 5]), f([6, 7])])
## whereas in Maple 18:  [f([2, 3]), f([4, 5]), f([6, 7])]
So it seems recursing is being done since 2015.
## Can it be prevented? Or is it a bug?
## I shall submit a bug report just in case ... .

Do I understand you correctly: You want to do some symbolic manipulations in Maple (maybe turn the system into a first order system) and then use Matlab to do the numerics?

If so you could use DEtoools[convertsys] in Maple to get a first order system. But why not start with doing the whole thing in Maple and let it worry about the conversion?

@nrebman1 Take a close look at the output from

diff(x[0, 1](t), t, t) = force(0)[1];

What did you expect on the right hand side?

There was an answer by Kitonum:

http://www.mapleprimes.com/questions/211685-Intersection-Of-Real-Intervals

@torabi Quite right. I overlooked that Q[4,2] is a sum:
200.*Y-2/(0.1e-1*y+.11)

thus won't be dealt with by py.
Hopefully this one works as intended:

px:=proc(u1) global X,x; local u:=expand(u1);
  if type(u,`+`) then return map(procname,u) end if;
  if not type(u,`*`) then return u end if;
  if member(X,{op(u)}) then diff(u/X,x) else u end if
end proc;
py:=subs(X=Y,x=y,eval(px));
pz:=subs(X=Z,x=z,eval(px));

@iman The problem is that all the highest derivatives are in the first equation.

You can solve that by doing:

sys:=select(has,dsys4,diff);
solve(sys,{diff(phi(x),x$3),diff(u(x),x$3),diff(w(x),x$4)}); #Cannot be done
##Now differentiate:
sys2:=diff(sys[2],x);
sys3:=diff(sys[3],x);
bcs23:=subs(x=0,convert(sys[2..3],D)); #New additional boundary conditions
sysnew:={sys[1],sys2,sys3}; #The new system
solve(sysnew, {diff(phi(x),x$3),diff(u(x),x$3),diff(w(x),x$4)}); #Now successful

@Jean-Michel Collard Although I probably won't have the need for using the CUDA capabilities (if any) of my computer I tried your example (i.e. the one from the help page) on my computer.
I got the same error message as you did.
At the same time a message came up from the graphics driver that it had stopped working.
I have NVIDIA GEFORCE GTX 860M version 364.72 on a 64 bit machine with Windows 10 and use Maple 2016.1 (64 bit).
Going to the NVIDIA page to look it seemed that it ought to have worked, but what do I know?
##
Just tried in Maple 2015.2 (64 bit) on the same machine: Same result.

I guess I don't need to tell you that you should avoid 2D math input as the plague!

Use 1D math input aka Maple input.

Well, if your ode system is 'secret' how do you imagine that we can help you?

This site is dedicated to questions about Maple.

@acer To check your result I tried an iteration approach treating the integral as a perturbation, but not replacing Int with anything else.
First letting f3(x) be zero in the integral, producing a result res[0]
Inserting the newly found f3 into the integral, compute again to find res[1], etc.
I stopped after 8 iterations with res[8].

The maximal differences between the results res[8] and your sol were:
9*10(-12), 5*10^(-13), and 5*10^(-9) for f1, f2, and f3, respectively.
The differences were larger for res[i] with i < 8, which seems to confirm that your result sol is pretty good.

For a very long time (since Maple 6) the new packages in Maple have been modules. LinearAlgebra was introduced in Maple 6 and as a module. Later some old ones have been turned into modules too.

Examples:
type(DEtools,`module`); #false
type(DEtools,table); # true
type(PDEtools,`module`); #true
type(plots,`module`); #true
Those packages were introduced before Maple 6 as tables.

@YasH Yes, automatic simplification is prevented in computing your a:

'5^( evalf(1.25, 30) )';
                                  

@torabi I suppose that I shall ignore the delta in front of all quantities?

I take the vertical bar |_x=0,1  to mean evaluated at x=0 and also at x=1, so there are two conditions for each equation.

Now the two second order conditions at x=0 and theta=0 (respectively) should be written

D[1,1](w)(0,theta)=0 and D[2,2](w)(x,0)=0.

D[1]    means the first order derivative w.r.t. the first variable (here x).
D[1,1] means the second order derivative w.r.t. the first variable (here x).
D[2]   means the first order derivative w.r.t. the second variable (here theta).

@marc sancandi Please notice that I revised the solution procedure. It had an unfortunate error. Neither of X[1] or X[2] should be applied to t. This mistake doesn't produce wrong results, but prevents the use of hardware floats (see below).
The computation with rosenbrock as well as with adamsfull now takes about 0.4s.

##
The following shows that the erroneous t in solproc prevents the use of hardware floats:

3(2); #Result 3. In Maple numbers can act as constant functions.
evalhf(3(2));
   
Error, unable to evaluate function `3` in evalhf