@Carl Love Yes, I wouldn't have thought of that, but isolating
 zip(solve, [(x-2)/3 , (y-1)/4 , (z-3)/3] =~ t, [x,y,z]);
made me understand.

@vv Certainly, yes. If this it not desirable then the first method should be chosen, i.e.

f:=proc(theta) proc(x) theta*x end proc end proc;

You have an initial value problem for the 1D wave equation. The solution is given by the result of
pdsolve({Eq1,u(x,0)=f(x),D[2](u)(x,0)=g(x)});

This is THE solution. So what more do you want?

@maple I edited my original answer above to include a shooting method. Take a look at that.
The conclusion I draw from that is that the original bvp solution is pretty good.

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);