The first parameter c0 in particular seems to give problems if it gets too large. Your comparison value is 5.7*10^(-6), but you use an interval c0=0..1.
Notice (incidentally) that y2(t) will be identically 0, so may be left out entirely.
For clarity I rewrote your procedure err (doesn't change the basic behavior).

Y1:=subs(sol, y1(t)): #This changes when parameters are changed!
err := proc (c0, n, s0, ks, h1, h2, kp, A, B) local sv1, sv2;
  sol(parameters = [c0, n, 175, s0, ks, h1, h2, kp, A, B]);
  sv1 := [Y1(1), Y1(100), Y1(210), Y1(2500), Y1(2800), Y1(3000)];
  sol(parameters = [5.7/10^6, 10.186, 175, 200, 1/20000000, 10000, .269, 1.5/10^7, 1.5, 2]);
  sv2 := [Y1(1), Y1(100), Y1(210), Y1(2500), Y1(2800), Y1(3000)];  
  add((sv1[i]-sv2[i])^2, i = 1 .. 6)
end proc;
err(5.7/10^6, 10.186, 200, 1/20000000, 10000, .269, 1.5/10^7, 1.5, 2);#Reference, 0.
err(2/10^4, 10, 200, 0, 10, .2, 1.5, 1.5, 2); #The value of ks is irrelevant
#I don't have access to GlobalOptimization, but tried the following, where the c0-interval is much smaller than yours:
Optimization:-NLPSolve(err,0 ..1e-4, 1 .. 20, 150 .. 250, 0 .. 1, 100 .. 15000, 0 .. .5, 0 .. 1, .5 .. 2, 1 .. 5);

I tried a shooting method: Finding the values for f''(0), g'(0), and theta'(0) that will result in f'(blt)=0, g(blt)=0, theta(blt)=0. I used fsolve for that.
I cheated by setting blt to 3, since it seemed from experimentation that f', g, and theta becomes zero rather rapidly.
Notice that I dont assign to blt directly.
restart;
with(plots):
k := .1; E := 1.0; Pr := 7.0; Ec := 1.0; p := 2.0;
Eq1 := diff(f(eta), eta, eta, eta)+f(eta)*(diff(f(eta), eta, eta))+Gr*theta(eta)-k*(diff(f(eta), eta))+2*E*g(eta) = 0;
Eq2 := diff(g(eta), eta, eta)+f(eta)*(diff(g(eta), eta))-k*g(eta)-2*E*(diff(f(eta), eta)) = 0;
Eq3 := diff(theta(eta), eta, eta)+Pr*(diff(theta(eta), eta))*f(eta)+Pr*Ec*((diff(f(eta), eta, eta))^2+(diff(g(eta), eta))^2) = 0;
bcs1 := f(0) = p, (D(f))(0) = 1, g(0) = 0, theta(0) = 1, theta(blt) = 0, (D(f))(blt) = 0, g(blt) = 0;
#The following gives a result for me, but Gr=9.5 does not.
res:=dsolve(eval({Eq1, Eq2, Eq3, bcs1}, {blt=3,Gr = 9.4}), numeric, output = listprocedure);
odeplot(res,[[eta,diff(f(eta),eta)],[eta,g(eta)],[eta,theta(eta)]],0..3);
#We take note of the results at eta=0. We need f''(0), g'(0), theta'(0) for later use:
res(0);
bcs1; #We shall use part of the info here in these:
ics := f(0) = p, D(f)(0) = 1, g(0) = 0, theta(0) = 1, D(theta)(0) = theta1, (D@D)(f)(0) =f2, D(g)(0) = g1;
respar:=dsolve({Eq1, Eq2, Eq3, ics}, numeric, output = listprocedure,parameters=[Gr,f2,g1,theta1]);
F1,G,TH:=op(subs(respar,[diff(f(eta),eta),g(eta),theta(eta)]));
#This outputs a list of results:
p:=proc(blt,Gr,f2,g1,theta1) option remember;
   respar(parameters=[Gr,f2,g1,theta1]);
   [F1(blt),G(blt),TH(blt)]
end proc;
blt0:=3: #Rather than 11.5
par:=[-.6,-1,-8]: #Taken from res(0)
#I inserted Gr0 = 11.5 to make it easier for fsolve to find a result.
#Note that fsolve wants a list of procedures, not a procedure outputting a list. Thus the need for p1,p2,p3.
for Gr0 in [10,11,11.5,12] do
  p1:=proc(f2,g1,theta1) p(blt0,Gr0,f2,g1,theta1)[1] end proc;
  p2:=proc(f2,g1,theta1) p(blt0,Gr0,f2,g1,theta1)[2] end proc;
  p3:=proc(f2,g1,theta1) p(blt0,Gr0,f2,g1,theta1)[3] end proc;
  par:=fsolve([p1,p2,p3],par);
  print(par);
  respar(parameters=[Gr0,op(par)]);
  q[Gr0]:=odeplot(respar,[[eta,diff(f(eta),eta)],[eta,g(eta)],[eta,theta(eta)]],0..blt0);
end do:
display(q[10],q[11],q[11.5],q[12]);

restart;
eq:=I*(a+3*b)+2*b+5*a = 3+2*I;
solve(evalc({Re,Im}~(eq)));

I thought I answered this very same question in
http://www.mapleprimes.com/questions/202715-Error-in-Dsolvenumericbvpconvertsys

The main problem was the range of r and theta. However, DEplot and dfieldplot don't require initial conditions. So to get something similar to the Matlab image you could do either

DEtools[DEplot](Sys1,[r(t),theta(t)],t=-1..1,r=0..7,theta=0..7);
## or
DEtools[dfieldplot](Sys1,[r(t),theta(t)],t=-1..1,r=0..7,theta=0..7);

Notice that the t-range is required, but not used.

Call your equation eq, then

eval(eq,t3=t1); #verifies that t3=t1 is a solution as you said.
params:=indets(eq,name) minus {t1,t3};
eq1:=simplify(eval(eq,params=~7)); #Setting all parameters to 7
res:=solve(eq1,t3); #Solving for t3
evalf[50](eval(res,t1=1));
#We see that t3 is significantly larger than t1 = 1.
Thus in general there is more one than one solution for t3. In special cases (e.g. the case params=~1) there is only one.
#You may try to find more than 2 solutions using fsolve
r2:=fsolve(eval(eq1,t1=1),t3,avoid={t3=1}); #Same as found from solve
fsolve(eval(eq1,t1=1),t3,avoid={t3=1,t3=r2},complex);
###Further explorations:
plots:-implicitplot(eq1,t1=-5..5,t3=-5..5,gridrefine=3);
plot(eval(lhs(eq1),t1=1),t3=-2..2,-1..1);
fsolve(eval(eq1,t1=1),t3=-1); #A negative real solution

Maybe try..catch could be useful:

printlevel:=2:
  for i from 0 to 3 do
    try
      1/i;
    catch:
    end try
  end do;

The following code works fine in Maple 18:
restart;
AIDS := [97, 206, 406, 700, 1289, 1654, 2576, 3392, 4922, 6343, 8359, 9968, 12990, 14397, 16604, 17124, 19585, 19707, 21392, 20846, 23690, 24610, 26228, 22768];
CAC := [seq(sum(AIDS[j]/(1000.0), j = 1 .. i), i = 1 .. 24)];
Time := [seq(1981+(i-1)*(1/2), i = 1 .. 24)];
#ln(CAC)=k*lnt+A
LnCAC := map(ln, CAC);
LnTime := map(ln, [seq((i+1)/(2*(1/10)), i = 1 .. 24)]);
with(stats);
fit[leastsquare[[x, y], y = k*x+lnA]]([LnTime, LnCAC]);
k := op(1, op(1, rhs(%))); LnA := op(2, rhs(%%)); A := exp(LnA);

#You may want to have a look at the new Statistics package. To make it do the exact same thing:
k:='k':
Statistics:-Fit(k*x+lnA,LnTime,LnCAC,x);
#Or use the original model directly:
Statistics:-Fit(t^k*A,Time,CAC,t);

The first problem is that from the 3 odes one cannot algebraically solve for the highest derivatives.
There is a way out of that:

restart;
#Your 3 odes:
sys:=[2*(diff(f1(x), x, x))+4*(diff(f2(x), x))+6*(diff(f3(x), x)) = 0, 10*f2(x)+12*(diff(f1(x), x))+14*f3(x) = 0, 16*(diff(f3(x), x, x, x, x))+19*(diff(f3(x), x, x))+22*(diff(f1(x), x))+25*f2(x)+56*f3(x)+63 = 0];
#Trying to solve for the highest derivatives (using solve):
solve(sys,{diff(f1(x),x,x),diff(f2(x),x),diff(f3(x),x$4)}); #No solution
#Solution: differentiate the second ode w.r.t. x:
ode:=diff(sys[2],x);
sys2:={sys[1],sys[3],ode};
#sys2 still has highest derivatives diff(f1(x),x,x),diff(f2(x),x),diff(f3(x),x$4).
#Checking solvability:
solve(sys2,{diff(f1(x),x,x),diff(f2(x),x),diff(f3(x),x$4)}); #A simple solution
#Your next problem is obvious: you have 12 boundary conditions. You can only have 2+1+4 = 7.
#You 5 conditions. I did the following.
bcs:={f1(1) = 0,f1(0) = 0 ,f2(1) = 0,f3(1) =0,f3(0) =0 , D(f3)(1) = 0,  D(f3)(0)=0};
res:=dsolve(sys2 union bcs,numeric,output=listprocedure);
plots:-odeplot(res,[x,f1(x)],0..1);

As pointed out in my earlier answer your initial conditions give problems. It turns out that you can get a solution on the closed interval [0, infinity) the following way:
restart;
fx:=piecewise(x<=0,0,x<1,1,0);
Eq := y(x)*diff(y(x), x, x)+diff(y(x), x)+y(x)=fx;
#Now isolate y''. Assume that y''(x)->q as x-> 0 from above.
#From y'' = (1-y-y')/y for x >0 (and x<1) we get by l'Hospital's rule using that we have an indeterminate form "0/0":
#q = limit((-y'-y'')/y',x=0,right) = -1-q. Thus q = -1/2.
#Now we may define Eqa as
Eqa:=simplify(diff(y(x),x,x)=piecewise(x=0,-1/2,(fx-y(x)-diff(y(x),x))/y(x))) assuming x>=0;
#We should expect that to be the valid replacement for Eq when x>=0.
res:= dsolve({Eqa, y(0) = 0, D(y)(0) = 1},numeric);
plots:-odeplot(res,[x,y(x)],0..5);
#There is no way that the solution can be extended to x<0 as a differentiable function. If you are satisfied with continuity at x=0 then y(x) = 0 is an extension.

pde:=diff(u(x,y),x,y)=sin(x+y);
bcs:=D[1](u)(x,0)=1,u(0,y)=(y-1)^2;
pdsolve({pde,bcs}); #No output
#We try mimicking what we would do by hand:
res:=pdsolve(pde); #General solution
eval(res,x=0);
eq1:=eval(%,{bcs});
eval(diff(res,x),y=0);
eq2:=eval(%,{bcs});
dsolve({eq2,eq1});
res2:=subs(%,res);

You got a problem with the initial condtions. Maple has to solve for the highest derivative i.e. y''. Thus you will be dividing by y. However, since one of your initial conditions is y(0)=0 you have a singularity at the initial point.
With y(0)<>0 there is no problem. I prefer using piecewise.

fx:=piecewise(x<=0,0,x<1,1,0);
Eq := y(x)*diff(y(x), x, x)+diff(y(x), x)+y(x)=fx;
ans := dsolve({Eq, y(0) = 1/2, D(y)(0) = 1},numeric);
plots:-odeplot(ans,[x,y(x)],-0.27..2);
#For these intial conditions there is a singularity at x = -0.2766....

If n is just a name, then
restart;
S:=sum(a[i]*x^i,i=0..n);
seq(S,n=0..4);
#Otherwise you could use mtaylor (instead of taylor to avoid the O-term):
restart;
n:=78;
S:=sum(a[i]*x^i,i=0..n);
mtaylor(S,x=0,7);#Notice the meaning of the last argument
seq(mtaylor(S,x=0,k),k=1..5);

Try using numerical integration.
Thus after having defined G and f (as functions) do:
w:=proc(r, phi, t) local xi,eta;
  if not type([_passed],list(realcons)) then return 'procname(_passed)' end if;
  evalf(Int(f(xi, eta)*G(r, phi, xi, eta, t)*xi, [xi = 0 .. 5, eta = 0 .. 2*Pi]))
end proc;
#You may now test w:
w(1, (1/5)*Pi, 0.1);
w(1, (1/5)*Pi, t); #Returning unevaluated by design
#Now plotting using the option adaptive= false to accomodate my limited patience:
plot(w(1, (1/5)*Pi, t), t = 0 .. .1, adaptive = false, numpoints = 25);

The behavior is the same in Maple 18 with simplex[convexhull] and geometry[convexhull].
You say it shouldn't be an issue of roundoff and yet it cannot be anything else. Both procedures convert the points to floats.