There is a syntax error in the bcs, but even after correcting it, pdsolve won't find the solution you gave.
It is, however, easy to check that the solution you give is in fact a solution.
 

PDE := diff(u(x, t), t) = (1/4)*exp(2-u(x, t))*diff(u(x, t), x, x)/(x^2+2);
ics := u(x, 0) = 2*(1-ln(2-x^2));
bcs := D[1](u)(0,t)=0, u(1,t)=2+ln(1+t); #syntax corrected
pdsolve({PDE,ics,bcs}); # No luck
## Checking your solution:
pdetest(u(x,t)=2 + ln(1+t) - 2*ln(2-x^2),[PDE,ics,bcs]);
## So your solution is indeed a solution to the problem.
### Now trying without conditions:
pdsolve(PDE); #getting two odes
sol:=pdsolve(PDE,build); #One ode solved, the other isn't.
ode_F1:=op([2,1],sol); #The unsolved ode
## We try dsolve:
dsolve(ode_F1); #No luck (and no surprise to me).
### Numerical solution. Notice slight change of syntax:
res:=pdsolve(PDE,{ics,bcs},numeric);
res:-animate(t=0..5); #An animation

## You can compare the animation produced above (frames=16) with an animation of the NASA solution:

plots:-animate(plot,[2 + ln(1+t) - 2*ln(2-x^2),x=0..1],t=0..5,frames=16);

So how was the solution produced in the NASA paper?

Without trying to understand why you are doing what you are doing, I can just observe that the integral in question is indeed divergent.
No way out of that.

Try in your worksheet:

Q1r:=convert(Q1,rational);
fct:=factor(Q1r*f);
eval(numer(fct),eta=1);#Not zero
denom(fct);

Since the denominator has the factor (1-eta) and the numerator is -5*A1*(-15*A1*A3+72*A2^2) at eta=1 (which quantity I shall assume is not zero) the integral int(Q1r*f, eta = 0 .. 1) is divergent.

Square brackets cannot be used as parentheses. Thus replace the [ ] in N with ( ).
Forget about the assumptions.
Use evalc(Re(...)) and evalc(Im(...)) instead of Re(...)  and Im(...).
evalc will see all names as real variables. So no need for assumptions as far as that is concerned.
On top you write Student(NumericalAnalysis). That has no effect.
Perhaps you meant with(Student:-NumericalAnalysis);

You must have set (inadventently) the output display to Maple Notation.
Change that by going to Tools/Options/Display/Output display and change to 2-D Math Notation.
Save globally before you leave.
 

This may not be elegant, but it works:
 

RN:=proc(yi::indexed) local ind1,ind2; 
  ind1:=op(1..-2,yi); 
  ind2:=op(-1,yi); 
  op(0,yi)[ind1,ind2+1] 
end proc;
RN(y[a,2]);
RN(y[a,b,c,2]);

Note. It could be written slightly simpler if you make use of the fact that the index is a sequence of only length 2.

I don't get any Maple errors when executing the code. I did notice, however a couple of syntactical problems, which were inconsequential though.
1.  It should be with(CurveFitting), not With(CurveFitting)., but you don't use that package anyway.
2. You write Digits; 5;  This has no effect at all.
3. It is a very bad idea to use x and x[i] in the same worksheet for different purposes. It can have bad consequences. Use two different names. Here it doesn't appear to be a problem though.
4. Use mul instead of product for concrete products.
6. Use add instead of sum for concrete sums.
###
The errors you are talking about, must be algorithmic.
I think you may want to change n-1 to n in a couple of places:
So your code without any changes to x and x[i] and printed more legibly and with a few corrections:

restart;
f := abs;
n := 8; 
h := 2/n; 
for i from 0 to n do x[i] := h*i-1 end do; 
for i from 0 to n do 
   L[i] := mul((x-x[j])/(x[i]-x[j]), j = 0 .. i-1)*mul((x-x[k])/(x[i]-x[k]), k = i+1 .. n) 
end do; 
lagr := add(f(x[p])*L[p], p = 0 .. n); 
P := expand(lagr); 
plot([P, f(x)], x -1 .. 1, y = 0 .. 2); 
M := maximize(abs(abs(x)-P), x = -1 .. 1); 
M1 := numapprox:-infnorm(abs(abs(x)-P),x=-1..1);

Try evalf(M);  after having executed the assignment to M.
Here are som plots. I show two of them together:

plot(P-f(x),x=-1..1,caption="The error");
plot([P, f(x)] ,x = -1 .. 1, y = 0 .. 2); p1:=%:
plot([seq([x[i],f(x[i])],i=0..n)],style=point,symbol=solidcircle, symbolsize=20,color=blue); pts:=%:
plots:-display(p1,pts);

If you do the integral in two steps: The inner first, then the outer, you get 3/sqrt(Pi).
For simplicity I do the conversion first, but that is not essential:

restart;
f:=x->1/sqrt(2*Pi)*exp(-x^2/2);
F:=convert(f(x)*f(y)*x*x*abs(x+y),piecewise,x);
Fy:=int(F,x=-infinity..infinity);
int(Fy,y=-infinity..infinity); #3/sqrt(Pi)
int(F,[x=-infinity..infinity,y=-infinity..infinity]); #Still incorrect

## I submitted an SCR.
## If you reverse the order of integration (which is clearly OK) then the iterated integral still gives the correct result, but the double range version does not, but returns 2/sqrt(Pi).
 

restart;
f:=x->1/sqrt(2*Pi)*exp(-x^2/2);
F:=convert(f(x)*f(y)*x*x*abs(x+y),piecewise,x);
Fx:=int(F,y=-infinity..infinity);
int(Fx,x=-infinity..infinity); # 3/sqrt(Pi)
int(F,[y=-infinity..infinity,x=-infinity..infinity]); # 2/sqrt(Pi)

The iterated integrals in both orders (and each in one command) works correctly too:
 

int(int(F,x=-infinity..infinity),y=-infinity..infinity);
int(int(F,y=-infinity..infinity),x=-infinity..infinity);

I was under the impression that the integral int(F, [x=-infinity..infinity,y=-infinity..infinity]) was just a convenient syntax and was handled as int(int(F,x=-infinity..infinity),y=-infinity..infinity).
Apparently I was mistaken.
In the help page for int, we find this note:

Note: The notation int(expression, [x = a..b, y = c..d]) is equivalent to int(int(expression, x = a..b), y = c..d) except that the single call to int accounts for the range of the outer variables (via assumptions) when computing the integration with respect to the inner variables.

So I guess that the handling of these assumptions went wrong somewhere.

2. To get the value of f'(eta) from your procedure p for, say, eta=2.3456, just do eval(diff(f(eta),eta),p(2.3456));
1. To plot results from different beta values, save the plots in say plt1,plt2,plt3, ... and then just use display from the plots package:
    display(plt1,plt2,plt3,...);

## If you need to find values for many values of eta, I suggest you use output=listprocedure in dsolve.
Then do this instead of the method above:
F1:=subs(p,diff(f(eta),eta));
F1(2.3456);
F1(0.76543);

I have several comments.
1. Why should different sets of equations have the same solutions?
  Consider the following very simple example, where the 3 equations were inspired by your l1,l2,l3:
 

eq:=exp(I*x)=1+I;
eqR:=evalc(Re(eq));
eqI:=evalc(Im(eq));
solve(eq,x,allsolutions);
evalc(%);
solve(eqR,x,allsolutions);
solve(eqI,x,allsolutions);

The solutions in this example are certainly not the same.

2. You need to use fsolve in all the 3 cases. You have solve in the last; it might be a misprint; you might have used fsolve in fact. I don't get any answer when using solve.
3. fsolve will only give you one answer, not all answers. Thus 3 different results might (in principle) all be correct.
 Check by using eval(eqs,result);
4. For your system(s) you should work with a high setting of Digits. I tried Digits:=30. You may want to try even higher.
5. When using evalc(expr) all variable names are assumed to represent real numbers.
6. Since you have to use fsolve, v[0] must be assigned a concrete value (as you also did).

I tried all 24 permutations and it seems that the order is consistently what you said:
The rectangles are plotted first and the curves second (i.e. on top). Furthermore the orders are opposite within the two groups.

 restart;
with(plots):with(plottools):
p1:=plot( x^3,x=-1..1,thickness=20,color=red):
p2:=plot(-x^3,x=-1..1,thickness=20,color=blue):
p3:=display(rectangle([0.5, 1],[0.75,-1],color=green)):
p4:=display(rectangle([-1,-0.1],[1,0.1],color=yellow)):
#####
L:=[p1,p2,p3,p4]:
LE:=ListTools:-Enumerate(L):
LE2:=subs({1=red,2=blue,3=green,4=yellow},LE):
LEP:=combinat:-permute(LE2,4):
for LL in LEP do 
  perm:=map2(op,1,LL);
  display(map2(op,2,LL),caption=typeset(perm))
end do;

I also tried the following to see the order of the arguments to the displayed plots:

disp:=display(L):
op(disp);
L2:=[p2,p4,p1,p3];
disp2:=display(L2);
op(disp2);

From this it appears that the order within L and disp is the same, and so is it for L2 and disp2.
So the answer probably cannot be found from examining Maple data structures, but must come from knowledge about the plot interface.
 

Try this where I have made life easy for myself in various ways.
 

restart;
plots:-animate(plot,[sin(a+x),x=-5..5],a=1..25); p:=%:
ppp:=indets(p,specfunc(CURVES)):
nops(ppp);
tr:=plottools:-transform((x,y)->[x,y,aa]); # aa unassigned
subs(aa=3,eval(tr)); #test to see if I succeed in planting 3 instead of aa in the procedure tr.
plots:-display(seq(subs(aa=i,eval(tr))(ppp[i]),i=1..25),insequence);

##Using existing lists
Repltlist3:=map(L->[op(L),t],Repltlist(t));
Impltlist3:=map(L->[op(L),t],Impltlist(t));
animate(pointplot3d,[Repltlist3,coords=cylindrical,symbol=solidsphere,symbolsize=20,color=red],t=-Pi..Pi,frames = 100,labels=[x,y,t]); P3:=%:
animate(pointplot3d,[Impltlist3,coords=cylindrical,symbol=solidsphere,symbolsize=20,color=blue],t=-Pi..Pi,frames = 100,labels=[x,y,t]); Q3:=%:
display(P3,Q3);

Your new ode is very closely connected to the one in your first question:
http://mapleprimes.com/questions/220659-Why-Maple-Gives-numeric-Exception

The ode in that question is:
ode1:=(a^2*x+(x^2-y(x)^2)*y(x))*diff(y(x),x)+x*(x^2-y(x)^2) = a^2*y(x);
## Your new ode is
ode2:=x*(a^2*x+(x^2-y(x)^2)*y(x))*diff(y(x),x)^2-(2*a^2*x*y(x)-(x^2-y(x)^2)^2)*diff(y(x),x)+a^2*y(x)^2-x*y(x)*(x^2-y(x)^2) = 0;
## If you solve for diff(y(x),x) as vv is doing in his answer:
odes:=solve(ode2,{diff(y(x),x)});
## then you get two odes one of which is the exact same as ode1 after that equation is solved for diff(y(x),x) too:
solve(ode1,{diff(y(x),x)});

It is therefore no surprise that you receive the same error in the two cases.

Obviously, this shouldn't happen.
Here is a workaround:
 

restart;
ode:=(a^2*x+(x^2-y(x)^2)*y(x))*diff(y(x),x)+x*(x^2-y(x)^2) = a^2*y(x);

ode_v:=PDEtools:-dchange({y(x)=sqrt(v(x))},ode,[v(x)]);
res:=dsolve(ode_v);
eval(res,v(x)=y(x)^2);
sol:=simplify(%) assuming y(x)>0;
odetest(sol,ode);

The result sol is given on implicit form as
sol := -x^2-2*a^2*arctanh(x/y(x))-y(x)^2+_C1 = 0

If you look for negative solutions instead, you get the same equation sol:
 

ode_v2:=PDEtools:-dchange({y(x)=-sqrt(v(x))},ode,[v(x)]);
res2:=dsolve(ode_v2);
eval(%,v(x)=y(x)^2);
simplify(%) assuming y(x)<0;

I shall submit an SCR on this.

Use length.

Examples:

M:=n->Matrix(n,symbol=m);
M(8);
length(M(8));
length(m[1,1]);
length(LinearAlgebra:-Determinant(M(9)));
length(LinearAlgebra:-Determinant(M(8)));
LinearAlgebra:-Determinant(M(9));
length(%);

As a response to the next but last command I get `[Length of output exceeds limit of 1000000]`, which is indeed the case according to the last command.