@Christian Wolinski I have no current access to Maple 5, but the behavior in Maple 8 is the same as it is in Maple 2018.

You really should use LinearAlgebra:-MatrixExponential(A,t) instead of DETools:-matrixDE .
The latter uses the very old matrix structure dating back to before Maple 6 (!).

@sand15 For A,B, and C try:
 

addressof(eval(A));
addressof(eval(B)); # As above for A
addressof(eval(C)); # Different

So C evaluates to a different copy of the table than A and B, which presumably is the purpose of 'copy'.

@Carl Love In the worksheet the contribution from phi disappears because the parameter a is set to zero in the loop.
The parameter a is taken from
 

Va := [0, 0];

The worksheet executes without any problem and as presented on my computer with Maple 2018.1.
## It works the same in Maple 2017.3 and also Maple 12.02.
## I'm using the 64 bit version on Windows 10.

@Carl Love A comment about results from dsolve/numeric/bvp spurred by your reply above.
By using the approximate solution
 

approxsoln=[f(eta)=lambda,diff(f(eta),eta)=eta/8-1,diff(f(eta),eta,eta)=0,g(eta)=0,theta(eta)=1]

I got results in the special case that were extremely good compared to the known exact solution in the special case K=0.
Using the same approximate solution and with abserr=1e-12 I got for Table 3, column 1:
 

[2.99473305260864, 3.07341413295920, 3.20904849825216, 3.38805933695604, 3.59894177889025]

The values reported in the paper are:

[2.98754,3.06345,3.19543,3.36981,3.57500]

Thus the differences are:
 

[0.7e-2, 0.10e-1, 0.14e-1, 0.18e-1, 0.24e-1]

So I agree with you that the numbers given in that column are not accurate to 4 decimal places.

## The approximate solution used I got from turning the ode for f into a system of first order in the special case. I noticed that that made dsolve give the desired solution without a user provided approxsoln. The initial profile generated by dsolve/numeric/bvp was in the special case:
 

[f(eta) = 3, f1(eta) = (1/8)*eta-1, f2(eta) = 0]

when %infinity=8, lambda=3 and the conversion was diff(f(eta),t)=f1(eta), diff(f1(eta),eta)=f2(eta).

@666 basha OK, so why not do that yourself?

@Carl Love Let us assume that F as presented indeed solves the ode for f with the given values of the parameters, in particular with lambda = 3. Then the given boundary condition f(0) = -lambda is WRONG and must be replaced by f(0) = lambda.
For that problem there are 2 solutions for f (at least 2). The first is found directly by dsolve/numeric, the second is found by dsolve/numeric with approxsoln=[f(eta)=F]:

restart:
Digits:= 15:
ODEs:= [
   #Eq 9/18/24:
   (1+K)*diff(f(eta),eta$3) + f(eta)*diff(f(eta),eta$2) - 2*n/(n+1)*diff(f(eta),eta)^2 -
   2/(n+1)*M*diff(f(eta),eta) + K*diff(g(eta),eta) + 2/(n+1)*sigma*theta(eta)  =  0,
   #Eq 10/19/25:
   (1+K/2)*diff(g(eta),eta$2) + f(eta)*diff(g(eta),eta) - 
   (3*n-1)/(n+1)*diff(f(eta),eta)*g(eta) - 2/(n+1)*K*(2*g(eta)+diff(f(eta),eta$2))  =  0,
   #Eq 11/20/26:
   diff(theta(eta),eta$2) + Pr*f(eta)*diff(theta(eta),eta)  =  0
]:   
params:={K= 0, Pr= 0.733, sigma= 0, lambda= 3, c= 1, M= 1, n= 1,%infinity=8};
odef:=eval(ODEs[1],params);
bcsf:=f(0) = lambda, D(f)(0) = -1, D(f)(%infinity)=0;
#Eq 32:
F:= simplify(eval[recurse](lambda - (1-exp(-eta*z))/z, [z= (lambda+sqrt(lambda^2+4*M-4))/2, lambda= 3, M= 1]));
odetest(f(eta)=F,{odef,bcsf}); # OK when lambda=3
res1:=dsolve(eval({odef, bcsf},params),numeric);
plots:-odeplot(res1,[eta,diff(f(eta),eta)]); p1:=%:
res2:=dsolve(eval({odef, bcsf},params),numeric,approxsoln=[f(eta)=F]);
plots:-odeplot(res2,[eta,diff(f(eta),eta)],color=blue); p2:=%:
plots:-odeplot(res2,[eta,diff(f(eta),eta)-diff(F,eta)]);
plots:-display(p1,p2);

The two solutions for diff(f(eta),eta):

The existence of two solutions could maybe be considered a consequence of the numerically necessary replacement of infinity by a finite number (here 8).
Notice that f(eta) = lambda - eta also satisfies odef and f(0) = lambda, D(f)(0) = -1, but not D(f)(infinity) = 0.
The acceptable solution of the two presented above must be the one for which f drops right away so that f ' (eta) appears as tending to zero as eta - > infinity.

@Carl Love There is a problem with the exact result in the special case where the ode for f can be solved separately.
In equation (32) in the paper I see a missing closing square bracket.
You interpret eq. 32 as:
F:= eval[recurse](lambda - (1-exp(-eta*z))/z, [z= (lambda+sqrt(lambda^2+4*M-4))/2, lambda= 3, M= 1]);
which seems reasonable to me. Thereby we get (after a trivial simplification):
F := 8/3+(1/3)*exp(-3*eta).
f(eta) = F solves the ode for f:

diff(f(eta), eta, eta, eta)+f(eta)*diff(f(eta), eta, eta)-diff(f(eta), eta)^2-diff(f(eta), eta) = 0

but not the first of the bcs since lambda = 3:

f(0) = -lambda, D(f)(0) = -1,D(f)(%infinity) = 0

So that is certainly a problem.
## I checked your code. I didn't find any difference between your version and the pdf file.

Since in Maple your example doesn't have a/b as a subexpression you may want to clarify your intention.
 

restart;
u:=y = s + (a/b)*log((a+c/b));
op([2,2],u);
op(%);
`*`(%);

The title you use I don't quite get. Is it important that the expression is an equation?

Yet another Dane with a corrupt file. Something is rotten in the state of Denmark.

@vv I just tried the conversion from 2D to 1D using lprint. It works. Thank you.

@AmirHosein Sadeghimanesh This is clearly a bug. Since I use 1D-input (aka Maple notation) , where double inequalities like 0 < x < 1 are not allowed I had to rewrite your input. I got the same results as you do though.
I tried without abs and got 29/360 with both piecewise and Heaviside.
I tried replacing abs(F) with sqrt(F^2) just in case, but (as they ought in this case) they gave the same result (although wrong).
 

restart;
q1:=-k5*x1*x2+k1*x1;
q2:=-k5*x1*x2+k2*x2;
p1:=piecewise(q1<=0,0,q1>=1,0,1);
p2:=piecewise(q2<=0,0,q2>=1,0,1);
F:=(k5*x2-k1)*(k5*x1-k2)-k5^2*x1*x2;
int(F*p1*p2,[k1 = 0 .. 1, k2 = 0 .. 1, k5 = 0 .. 1, x1 = 0 .. 1, x2 = 0 .. 1]); # 29/360
int(sqrt(F^2)*p1*p2,[k1 = 0 .. 1, k2 = 0 .. 1, k5 = 0 .. 1, x1 = 0 .. 1, x2 = 0 .. 1]); # infinity
int(abs(F)*p1*p2,[k1 = 0 .. 1, k2 = 0 .. 1, k5 = 0 .. 1, x1 = 0 .. 1, x2 = 0 .. 1]); # infinity
g:=convert(p1*p2,Heaviside);
int(F*g,[k1 = 0 .. 1, k2 = 0 .. 1, k5 = 0 .. 1, x1 = 0 .. 1, x2 = 0 .. 1]); # 29/360
int(sqrt(F^2)*g,[k1 = 0 .. 1, k2 = 0 .. 1, k5 = 0 .. 1, x1 = 0 .. 1, x2 = 0 .. 1]); # -29/360
int(abs(F)*g,[k1 = 0 .. 1, k2 = 0 .. 1, k5 = 0 .. 1, x1 = 0 .. 1, x2 = 0 .. 1]); #  -29/360

@Mariusz Iwaniuk It would be nice if your title in comments contained more than just a period. In the notifications from MaplePrimes the link to those comments are basically invisible. Until now I have thought there were none.

@Mariusz Iwaniuk There is a bug in convert/Heaviside, which affects your worksheet:
 

restart;
p:=piecewise(t <> 8*n, y(t), t = 8*n, 0);
pH:=convert(p, Heaviside);

The output from the conversion is   -y(t)*Dirac(t-8*n)+y(t), which is clearly not the same as p.
Just do:

int(eval(p,{y(t)= t, n=1}),t=0..10);
int(eval(pH,{y(t)= t, n=1}),t=0..10);

Results 50 and 42.

@mmcdara It appears that the bug in dsolve is due to a bug in convert/Heaviside:
 

restart;
ode:=diff(x(t),t)=piecewise(t=1,0,x(t));
convert(ode,Heaviside);

Result: diff(x(t), t) = -x(t)*Dirac(t-1)+x(t), which is clearly not the same as ode.
This bug is at least as old as Maple 8.