When I run your worksheet in Maple 2015 without trying to understand what is going on I don't get any unevaluated integrals.
I see that you have Maple 13 and from the output you show above you have some instances of Int (grey integrals), i.e. the inert version of int.
Try value on Vc[1](t) to turn Int into int.
Maple 2015 gives you:

@maple What do you mean by a good method? Is there anything wrong with the dsolve/numeric result I gave above:

ode:=diff(f(x),x,x,x)+f(x)*diff(f(x),x,x)+1-diff(f(x),x)^2=0;
bcs:=f(0)=0,D(f)(0)=0,D(f)(infinity)=1;
res:=dsolve(eval({ode,bcs},infinity=5),numeric); #Change to infinity=whatever if you like.

@maple Try setting (as Tom Leslie does) infinity=30 in dsolve. Then plot the result only on the interval x=0..5. Save the plot in p30, say.
Then set infinity=5 in dsolve. Plot the result on the interval x=0..5. Save the result in p5.
Now do
plots:-display(p30,p5);

You should see that the one is on top of the other, i.e. there is no visible difference, and visible difference is what matters if you are plotting (obviously).

@acer You are right about the need for protecting possible calls to Int within u.
But Diff and diff do accept the variables in a list, in fact in at least one case it is needed:
diff( x^2, []);

A revised version:
CollapseNestDiff:=proc(u::specfunc({diff,Diff,%diff})) local a,b,INT;
   a:=subs(Int=INT,diff=Int,Diff=Int,%diff=Int,u);
   b:=IntegrationTools:-CollapseNested(a);
   subs(Int=Diff,INT=Int,b)
end proc;
#One example:
u:=Diff(Diff(Diff(x^2*exp(x*y)*Int(sin(y)*x^3,x),x),x),y);
lprint(u);
value(%);
CollapseNestDiff(u);
lprint(%);
value(%);

@acer In the meantime, if there is a strong enough need for such a facility, then one could use IntegrationTools:-CollapseNested for this purpose:

CollapseNestDiff:=proc(u::specfunc({diff,Diff,%diff})) local a,b;
   a:=subs(diff=Int,Diff=Int,%diff=Int,u);
   b:=IntegrationTools:-CollapseNested(a);
   subs(Int=Diff,b)
end proc;

@John Fredsted That switching surprised me, but take a look at the lprint versions:

expr := Diff(f(x,y),x,y,y):
lprint(expr);
Diff(f(x, y), x, y, y)
expr1 := convert(expr,diff):
lprint(expr1);
diff(diff(diff(f(x, y), x), y), y)
expr2 := convert(expr1,Diff):
lprint(expr2);
Diff(Diff(Diff(f(x, y), x), y), y)

@Honigmelone You may also want to look at Physics:-diff:

restart;
a:=b(t)+c(t)*diff(b(t),t)+diff(b(t),t$2)+diff(b(t),t$3)+45*diff(c(t),t);
Physics:-diff(a,diff(b(t),t));
Physics:-diff(a,diff(c(t),t));
Physics:-diff(a,diff(b(t),t,t));

@iman What kind of symmetry are you thinking about? A graph of a function of x couldn't possibly be symmetric about the x-axis unless the function is zero on the whole interval.

The solution you are getting I wouldn't trust though.
Try doing:
SYS:=expand(solve(newsys,{diff(g1(y),y$3),diff(g2(y),y$4),diff(g3(y),y$4)}));
and observe the denominators on the right hand sides.

You can expect problems for y=0.
You set abserr=1e10, i.e. 10^10. That is kind of large for an absolute error where the values of g1,g2,g3 are in fact pretty small (~10^(-9)).

@iman L is never given a value. What is it?

@iman No in my two uses of dsolve above a finite difference method is used.
Shooting as I understand it works like this on a boundary value problem.
Again an example:
restart;
sys:={diff(y(x),x,x)+lambda(x)*y(x)=0,diff(lambda(x),x)};
bcs:={y(0)=0,y(Pi)=0,D(y)(0)=1};
##The finite difference result first:
res1:=dsolve(sys union bcs,numeric);
# Now a shooting method:
ics:={y(0)=0,D(y)(0)=1,lambda(0)=a}; #A parameter a to be determined
##Notice that ics are initial conditions, not boundary conditions.
## We must determine the parameter a such that y(Pi)=0:
res2:=dsolve(sys union ics,numeric, parameters=[a],output=listprocedure);
Y:=subs(res2,y(x)); #Procedure to evaluate y(x)
Q:=proc(a) res2(parameters=[a]); Y(Pi) end proc;
L:=fsolve(Q,0.1234); #even without the guess it works: fsolve(Q)
res2(parameters=[L]);
res2(0);
res1(0);

@iman
I believe that the method used by dsolve/numeric/bvp to find a parameter is the way I do it in the following simple example (the second way below). I know of no name for that method, but it is pretty straightforward.

restart;
ode:=diff(y(x),x,x)+lambda*y(x)=0;
res1:=dsolve({ode,y(0)=0,y(Pi)=0,D(y)(0)=1},numeric);
res1(0);
## Now the explicit way of doing the same:
sys:={subs(lambda=lambda(x),ode),diff(lambda(x),x)=0};
bcs:={y(0)=0,y(Pi)=0,D(y)(0)=1};
res2:=dsolve(sys union bcs,numeric);
res2(0);

@iman The method used used by dsolve/numeric/bvp is a finite difference method, not a shooting method.

Besides the problems mentioned by acer, it seems that no value is ever given to a, which appears in extra_bcs.

Switching to 1D input (as suggested by acer) will make it easier for me (and I believe you) to understand what is going on.

What was the intended meaning of this last  one of your set ODE:

diff(y(eta), eta)+Sb*f(eta)*y(eta)-Pe*(w(eta), eta)*(diff(t(eta), eta))-Pe*t(eta)*y(eta) = 0;

Instead of (w(eta), eta) did you want just w(eta) or did you want diff(w(eta), eta) ?

If I assume the former (i.e. w(eta) ) then by doing

etainf:=15:
BCS2:=BC union remove(hastype,IC,name);
dsolve(ODE union BCS2,numeric);
# I get the error message:
Error, (in dsolve/numeric/bvp/convertsys) too many boundary conditions: expected 12, got 15

@tomleslie I think the comment may have been meant for me, but put in a wrong place.
I tried to answer his question as he put it: "Problem with seq", i.e. as a problem having to do with seq.
Then he had the (to me) puzzling comment to my answer saying: "However, the output gives commas and quotes".
I now see from this latest comment that he is (now at least) using print and not printf.

@das1404 I have added an example below. Your comment "However, the output gives commas and quotes" puzzles me.