@satya123999 I tried your file with ImportMatrix and ImportVector:

A:=ImportVector("G:/MapleDiverse/MaplePrimes/test2.txt",datatype=string);
parse(A[1]);

Worked fine.
I first tried a text file test.txt containing several lines:

1.0*x1+1.0*x2+1.0*x3+1.0*x4+1.0*x5+1.0*x6+1.0*x7+1.0*x8+1.0*x9+1.0*x10+1.0*x11+1.0*x12
a+7*x+.78*y
45*x^2+sin(x)
Pi*u+8.98765+z

For that file I didn't need datatype=string:
V:=ImportVector("G:/MapleDiverse/MaplePrimes/test.txt"); #ImportMatrix works too.
V[1]; #outputs the algebraic expression in the first line above.

Finally I tried a text file like test2.txt, but with the difference that before I saved it I changed the line.
Now I didn't need datatype=string, but the output from
V:=ImportVector("G:/MapleDiverse/MaplePrimes/test3.txt");
was a vector with two elements V[2] being NULL.

restart;
A:=Int(x*cos(x)/(1+3*sin(n*x)^2), x = 0 .. Pi);
#Will split the integral in integrals over a..b with a=(2*j-1)/n*Pi/2, b=a+Pi/n, j=1..n-1.
#Since we are interested in limits we need not worry about the missing intervals at the beginning and the ends: 0..Pi/2/n and Pi-Pi/n .. Pi:
A1:=Int(x*cos(x)/(1+3*sin(n*x)^2), x=a .. b );
IntegrationTools:-Parts(A1,x*cos(x));
#We need to take limits from above a and below b, but this is more conveniently done like this:
A2:=subs(arctan(2*tan(n*b))=Pi/2,arctan(2*tan(n*a))=-Pi/2,%);
Bint,Btrig:=selectremove(has,A2,Int);
#The sum of the contributions from the integrated terms:
Sum(subs(b=a+Pi/n,a=(2*j-1)/n*Pi/2,Btrig),j=1..n-1);
S:=value(%);
limit(S,n=infinity);
#Now we need only show that the sum of the integral parts converges to zero:
#The individual integrals are each represented by
Bint;
#Their sum can be be combined, which we can do simply by:
B:=eval(Bint,{a=0,b=Pi});
#B is bounded in absolute value by
Int((1/2)*Pi/2*(1+x)/n, x = 0 .. Pi);
value(%);
#clearly tending to zero as n->infinity.

The best I could come up with was:

f:=s->Int(t^2/(1+s*cosh(2*t))^2,t=-infinity..infinity);
value(f(1)); #SOo s=1 works OK
series(f(s),s=1,10);

Dipping into the problem slightly deeper it seems that the problem is that the second derivative of H appears, but it can be eliminated by using eq7.
At the same time one might as well use eq8 to elimate h' from the other equations.

sys16a:=[seq(cat(eq,i),i=1..6)];
indets(sys16a);
H'' appears, but can be eliminated
sys78:=solve({eq7,eq8},{diff(H(eta), eta),diff(h(eta), eta)});
sys16b:=subs(sys78,sys16a);
indets(sys16b);
#Now you can solve for the highest derivatives (H'' is gone):
sys16:=solve(sys16b,{diff(F(eta), eta, eta), diff(G(eta), eta, eta), diff(P(eta), eta), diff(Q(eta), eta), diff(f(eta), eta), diff(g(eta), eta)});
#So the total system (excluding boundary conditions:
sys:={op(sys16)} union sys78;

You will still have problems with the boundary value problem.
However, to convince oneself that things are not all bad and worth pursuing further one might try (as I did) an initial value problem at eta=8:
ic:={ F(8) = 0, G(8) = 0, D(F)(8)=-0.1, D(G)(8)=.1,h(8)=0.1,H(8) = 1,f(8) = 0, g(8) = 0, Q(8) = .2, P(8) = 0};
res:=dsolve(eval(sys,para) union ic,numeric);
L:=[[eta,F(eta)],[eta,G(eta)],[eta,g(eta)],[eta,f(eta)],[eta,Q(eta)],[eta,P(eta)]];
for i to 6 do p[i]:=plots:-odeplot(res,L[i],5.1..8) end do:
plots:-display(Array([[p[1],p[2],p[3]],[p[4],p[5],p[6]]]));

convertsys needs to solve for the highest derivatives of the unknown functions. It simply cannot do that.

Highest_deriv:={diff(F(eta), eta, eta), diff(G(eta), eta, eta), diff(H(eta), eta, eta), diff(P(eta), eta), diff(Q(eta), eta), diff(f(eta), eta), diff(g(eta), eta), diff(h(eta), eta)};

solve(subs(para, {eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8}),Highest_deriv);

returns NULL, meaning no solution.
Thus you need to think a little more about your system.

It really is just an algebraic problem. Consider all functions as just names. Tru to solve for the highest derivatives:

HD:=freeze~(Highest_deriv); #Freezing the highest derivatives: The unknowns
res1:=subs(Highest_deriv=~HD,{eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8});
fu:=freeze~(indets(res1,function)); #Freezing the other functions
res2:=subs(indets(res1,function)=~fu,res1); #System in names only
solve(res2,HD); #No solution

The matrix is 6x6, so the characteristic polynomial is of order 6. Thus only in special cases can you expect to find its roots explicitly because of the variable test10.

Example:

close3:=x->evalf(x); #Trivial example for close3
t:=.12345; #Using also a concrete value for t.
NewMatrix3 := Matrix([[test10, close3(t), close3(t+1), close3(t+2), close3(t+3), close3(t+4)], [close3(t), close3(t+1), close3(t+2), close3(t+3), close3(t+4), close3(t+5)], [close3(t+1), close3(t+2), close3(t+3), close3(t+4), close3(t+5), 0], [close3(t+2), close3(t+3), close3(t+4), close3(t+5), 0, 0], [close3(t+3), close3(t+4), close3(t+5), 0, 0, 0], [close3(t+4), close3(t+5), 0, 0, 0, 0], [close3(t+5), 0, 0, 0, 0, 0]]);

A:=LinearAlgebra:-Transpose(NewMatrix3).NewMatrix3;
LinearAlgebra:-Eigenvalues(A);
No reason to attempt finding eigenvectors when you are already stuck here.
However, for concrete values of test10 there is no problem:

LinearAlgebra:-Eigenvalues(eval(A,test10=7.913));
simplify(%);
LinearAlgebra:-Eigenvectors(eval(A,test10=7.913));
simplify([%]);

Here is a solution done in part like by hand:

k := diff(u(t, x), t,t) = diff(u(t, x), x,x);
bc := u(0, x) = 1/(1+x^2);
v:=D[1](u)(0,x)=0; #You need to use D, not diff.
res:=pdsolve(k);
eq1:=subs(t=0,bc,res);
diff(res,t);
convert(%,D);
eq2:=subs(t=0,v,%);
eq3:=map(int,eq2,x)+(C=0);
solve({eq1,eq3},{_F1(x),_F2(x)});
_F1,_F2:=op(unapply~(rhs~(%),x));
res;
sol:=normal(%);

What first comes to my mind is mtaylor:

pol:=x^5 + 4*x*y^4 + 2*y^3 +  x*y^2 + x^2 + y + 3;
mtaylor(pol,[x,y],4);

Another way:
map(select,q->evalb(degree(q)<=3),pol);

You may want to take a look at the CurveFitting package.

Example:
restart;
X:=[3,4,5,6,7];
Y:=[7.42494922444550, 3.67768248674133, 2.52235142453921, 1.95610223891559, 1.61770309810016];
f:=CurveFitting:-Spline(X,Y,x);
ptplot:=plot(X,Y,style=point,symbolsize=20,color=blue,symbol=solidcircle);
plot(f,x=3..7); cplot:=%:
plots:-display(cplot,ptplot);

Extrapolation is a dangerous thing to do, but if you have a model function you could try using CurveFitting:-LeastSquares as I now see that Kitonum is doing (you had two versions of the question, I just deleted the one to which I first answered).

I don't think there is any such thin.
Here is an ad hoc solution:

eps:=.05:
plot3d(`if`(abs(y^2+cos(x))<eps,undefined,(x^2+sin(x))/(y^2+cos(x))), x = -7 .. 7, y = -5 .. 5, view = [-7 .. 7, -6 .. 6, -1 .. 7], numpoints = 50000, scaling = constrained, axes = normal);

A variation of the same theme:

u:=(x^2+sin(x))/(y^2+cos(x)):
plot3d(piecewise(abs(u)>100,undefined,u), x = -7 .. 7, y = -5 .. 5, view = [-7 .. 7, -6 .. 6, -1..7], numpoints = 50000, scaling = constrained, axes = normal);

A print statement produces no output (i.e. NULL). The last statement in a procedure determines the output (unless there is an explicit return statement somewhere).
When a call to a procedure is made the formal parameters are just replaced with their values. No more evaluation.
Thus your procedure generate_x when called as in
generate_x(m);
m;
#or
generate_x(f(t));
f(t);
will return an integer between 1 and 10 and work as long as m and f(t) are assignable objects.
#########
Your procedure generate_y produces no output (NULL) since the last statement is print(y).
Thus
generate_y(k);
#prints k (literally), but produces no output. Try
%; #Shows the last output: the procedure definition (if that was executed just before the call to generate_y).
But the call did assign to k:
k;
Replacing print(y) with NULL has the same effect, but without printing of course.

I just tried in Maple 15. I got 'pt' just like you. Apparently a bug of some sort. pt is a local variable.
The quartic polynomial has 4 roots of course. I don't know if that is part of the problem.
But try this:
restart;
A := solve(x^4 + x^3 + x^2 + x + 1  = n, x,explicit ): #Notice 'explicit'
evalf(eval([A],n=100000));
asympt( A[3], n, 2 );
asympt( A[4], n, 2 );


In Maple 17.02 there is no such problem:

restart;
A := solve( x^4 + x^3 + x^2 + x + 1  = n, x );
asympt( A, n, 2 );
                            

asympt( A, n, 3 );
                            

Which version are you using?

To test the validity of the expansion I continued with:
res3:=convert(%,polynom);
N:=10000:
fsolve( x^4 + x^3 + x^2 + x + 1  = N, x );
evalf(eval(res3,n=N));
#Reasonable agreeness.

There is a problem to begin with. The differential equation is a cubic equation in the leading derivative (in this case simply the first).
Thus using the usual approach dsolve will have to solve for that derivative. A cubic has 3 roots. They may or may not all be real. 
In your case there is one real root just assuming real parameters and beta>0 suffices.

restart;
diffeq := diff(w(r), r)+2*beta*diff(w(r),r)^3-(1/2)*S*(r-m^2/r) = 0;
con := w(1) = 1;
pol:=subs(diff(w(r),r)=w1,lhs(diffeq)); #The polynomial in the derivative w1
d1:=discrim(pol,w1); #Its discriminant
d:=collect(d1,beta,factor);#We se that beta>0 suffices for one real root since then d<0.
#Your parameters:
param:={beta=0.2,S=0.5,m=2};
#Illustration of the comments above:
fsolve(eval(pol,param union {r=3}),w1); #One real root
fsolve(eval(pol,{beta=-.01,S=1,m=1,r=3}),w1); #3 real roots
#We solve for the highest derivative:
res:=solve(diffeq,{diff(w(r),r)}): #3 solutions
res:=simplify([res]) assuming beta>0; #Simplifying
eval(res,param union {r=3}); #Testing to see which is real
evalf(%);
#Now trying res[1]
sol:=dsolve(eval({op(res[1]),con},param),numeric);
plots:-odeplot(sol,[r,w(r)],0.001..10);
#Seems OK.
###############
## The following is not so straightforward, but has definite advantages. For help see
?dsolve,numeric,ivp
#First we take your concrete polynomial:
pol1:=eval(pol,param);
#The we define a procedure to be fed to dsolve. This procedure solves the cubic numerically using fsolve, which only returns real solutions:
p:=proc(N,rr,w,wp)
   wp[1]:=fsolve(eval(pol1,r=rr),w1);
end proc;
sol2:=dsolve(numeric,number=1,procedure=p,start=1,initial=Array([1]),procvars=[w(r)]);
plots:-odeplot(sol2,[r,w(r)],0..10);
#Looks good!
#Now a version which should work for all parameters for which pol has precisely one real root.
p2:=proc(N,rr,w,wp,beta1,S1,m1)
   wp[1]:=fsolve(eval(pol,{r=rr,beta=beta1,S=S1,m=m1}),w1);
end proc;
#Notice that the extra parameters to p2 are the values of beta,S, and m in that order.
sol3:=dsolve(numeric,number=1,procedure=rcurry(p2,0.2,0.5,2),start=1,initial=Array([1]),procvars=[w(r)]);
plots:-odeplot(sol3,[r,w(r)],0..10);
#Looks good too!

In a vain attempt to speed up things I tried turning the bvp into an initial value problem with parameters. That would involve fsolve so the computation overall took about twice as much time.

However, in the proces I discovered that the boundary value problem can have at least two solutions for given values of a,b,c,d.This is certainly the case for a=b=c=d=0.1 as the following shows. The two solutions have different values for f''(0).

restart;
eq1:= diff(f(eta),eta$3) - a*diff(f(eta),eta$1)^2 + b*f(eta)*diff(f(eta),eta$2) = 0:
bc:= f(0)=0, D(f)(0) = 1+c*(D@@2)(f)(0), D(f)(8)=d:
sys:= unapply({eq1,bc}, a, b, c, d):
Sol:= (a,b,c,d)-> dsolve(sys(a,b,c,d), numeric):
#Example a = b = c = d = 0.1.
plots:-odeplot(Sol(.1,.1,.1,.1),[eta,f(eta)],0..8,caption="Solution 1"); p1a:=%:
#############
ic:= f(0)=0, D(f)(0) = 1+c*f2, (D@@2)(f)(0)=f2:
solpar:=dsolve({eq1,ic}, numeric,parameters=[a,b,c,f2],output=listprocedure):
F1:=subs(solpar,diff(f(eta),eta));
#We make a procedure returning f'(8)-d:
p:=proc(a,b,c,d,f2)
   if not type([_passed],list(numeric)) then return 'procname(_passed)' end if;
   solpar(parameters=[a,b,c,f2]);
   F1(8)-d
end proc;
plot(p(.1,.1,.1,.1,f2),f2=-1..1); #Two zeros!!
r1:=fsolve(p(.1,.1,.1,.1,f2)=0,f2=-.3);
r2:=fsolve(p(.1,.1,.1,.1,f2)=0,f2=-.8);
solpar(parameters=[.1,.1,.1,r1]):
plots:-odeplot(solpar,[eta,f(eta)],0..8,color=red,caption="Solution 1 again"); p1b:=%:
solpar(parameters=[.1,.1,.1,r2]):
plots:-odeplot(solpar,[eta,f(eta)],0..8,color=blue,caption="Solution 2"); p2b:=%:
#The second solution can also be produced by solving the bvp problem directly when an approximate solution is known, and of course we have one.

Sol2:= dsolve(sys(.1,.1,.1,.1), numeric,approxsoln=solpar);

plots:-odeplot(Sol2,[eta,f(eta)],0..8,caption="Solution 2 again"); p2a:=%:
plots:-display(p1a,p2a);
plots:-display(p1b,p2b);
plots:-display(p1a,p2a,p1b,p2b);

I don't quite understand what you mean by plotting epsilon against time.
But if the intention is to plot the dependence on epsilon of the solution y then you could do something like this:

restart;
dsys1 :=diff(y(t),t)=y(t)*((1-y(t)/3-epsilon)-0.8*y(t)/(y(t)^2+0.5^2));
ini1:=y(0)=0.5:
dsol1 :=dsolve({dsys1,ini1},numeric, output=listprocedure,range=0..1,parameters=[epsilon]);
dsolu:=subs(dsol1,y(t));
#Setting the parameter epsilon for a test plot:
dsol1(parameters=[.02]);
plots:-odeplot(dsolu,[t,y(t)],0..1,caption=typeset(epsilon = 0.02));
#Now set up a procedure which returns unevaluated for non-numeric input, but otherwise sets the parameter epsilon to the input value (eps) and simply returns the procedure dsolu.
p:=proc(eps)
  if not type(eps,numeric) then return 'procname(_passed)' end if;
  dsol1(parameters=[eps]);
  dsolu
end proc;
#An animation:
plots:-animate(plot,[p(epsilon)(t),t=0..1],epsilon=0..1);
#A 3d plot:
plot3d(p(epsilon)(t),t=0..1,epsilon=0..1);