The ODE system has two parameters. One, A, will get a fixed value. The
other, U0, will be used like an independent variable for plotting.
> 
eq1:= diff(r[1](t), t, t)+(.3293064114+209.6419478*U0)*(diff(r[1](t), t))
+569.4324330*r[1](t)0.3123434112e2*V(t) = 1.547206836*U0^2*q(t):
eq2:= 2.03*10^(8)*(diff(V(t), t))+4.065040650*10^(11)*V(t)
+0.3123434112e2*(diff(r[1](t), t)) = 0:
eq3 := diff(q(t), t, t)+1047.197552*U0*(q(t)^21)*(diff(q(t), t))
+1.096622713*10^6*U0^2*q(t) = 2822.855019*A*(diff(r[1](t), t, t)):

> 
ics:=r[1](0)=0,D(r[1])(0)=1e1,V(0)=0,q(0)=0,D(q)(0)=0:

> 
res := dsolve({eq1,eq2,eq3, ics},numeric,output=listprocedure,parameters=[A,U0]):

I will call the procedure returned by dsolve, for evalutions of V(t), as the
dsolve numeric solutionprocedure in the discussion below.
> 
WV(parameters=[A=1e0]):

The goal is to produce a 3D plot of V(t) as a function of t and the parameter U0.
> 
tlo,thi := 0.0, 2.0;
U0lo,U0hi := 1e3, 2e1;

This is the grid size used for plot3d below. It is nothing special.
First, I'll demonstrate that a 3D plot can be produced quickly, by populating a
Matrix for floatingpoint evaluations of V(t), depending on t in the first
Matrixdimension and on parameter U0 in the second Matrixdimansion.
The surfdata command is used. This is similar to how plot3d works.
This computes reasonably quickly.
But generating the numeric values for U0 and t , based on the i,j positions
in the Matrix, involves the kind of sequence generation formulas that are
error prone for people.
> 
str := time[real]():
M:=Matrix(m,n,datatype=float[8]):
for j from 1 to n do
u0 := U0lo+(j1)*(U0hiU0lo)/(n1);
WV(parameters=[U0=u0]);
for i from 1 to m do
T := tlo+(i1)*(thitlo)/(m1);
try
M[i,j] := WV(T);
catch:
# mostly maxfun complaint for t above some value.
M[i,j] := Float(undefined);
end try;
end do:
end do:
plots:surfdata(M, tlo..thi, U0lo..U0hi,
labels=["t","U0","V(t)"]);
(time[real]()str)*'seconds';

So let's try it using the plot3d command directly. A 2parameter procedure
is constructed, to pass to plot3d. It's not too complicated. This procedure
will uses one of its numeric arguments to set the ODE's U0 parameter's
value for the dsolve numeric solutionprocedure, and then pass along
the other numeric argument as a t value.
It's much slower than the surfdata call above..
> 
VofU0 := proc(T,u0)
WV(parameters=[U0=u0]);
WV(T);
end proc:

> 
str := time[real]():
plot3d(VofU0, tlo..thi, U0lo..U0hi,
grid=[m,n], labels=["t","U0","V(t)"]);
(time[real]()str)*'seconds';

One reason why the previous attempt is slow is that the plot3d command
is changing values for U0 in its outer loop, and changing values of t in its
inner loop. The consequence is that the value for U0 changes for every
single evaluation of the plotted procedure. This makes the dsolve numeric
solutionprocedure work harder, by losing/discarding prior numeric
solution details.
The simple 3D plot below demonstrates that the plot3d command chooses
xy pairs by letting its second supplied independent variable be the one
that changes in its outer loop. Each time the value for y changes the counter
goes up by one.
> 
glob:=0:
plot3d(proc(x,y) global glob; glob:=glob+1; end proc,
0..3, 0..7, grid=[3,3],
shading=zhue, labels=["x","y","glob"]);

So now let's try and be clever and call the plot3d command with the two
independent variables reversed in position (in the call). That will make
the outer loop change t instead of the ODE parameter U0.
We can use the transform command to swap the two indepenent
axes in the plot, if we prefer the axes roles switched. Or we could use the
parametric calling sequence of plot3d for the same effect.
The problem is that this is still much slower!
> 
VofU0rev := proc(u0,T)
WV(parameters=[U0=u0]);
WV(T);
end proc:

> 
str := time[real]():
Prev:=plot3d(VofU0rev, U0lo..U0hi, tlo..thi,
grid=[n,m], labels=["U0","t","V(t)"]):
(time[real]()str)*'seconds';
plots:display(
plottools:transform((x,y,z)>[y,x,z])(Prev),
labels=["t","U0","V(t)"],
orientation=[50,70,0]);

There is something else to adjust, to get the quick timing while using
the plot3d command here.
It turns out that setting the parameter's numeric value in the
dsolve numeric solutionprocedure causes the loss of previous details
of the numeric solving, even if the parameter's value is the same.
So calling the dsolve numeric solutionprocedure to set the parameter
value must be avoided, in the case that the new value is the same as
the old value.
One way to do that is to have the plotted procedure first call the
dsolve numeric solutionprocedure to query the current parameter
value, so as to not reset the value if it is not changed. Another way
is to use a local of an appliable module to store the running value
of the parameter, and check against that. I choose the second way.
And plot3d must still be called with the first independent variablerange
as denoting the ODE's parameter (instead of the ODE's independent
variable).
And the resulting plot is fast once more.
> 
VofU0module := module()
local ModuleApply, paramloc;
ModuleApply := proc(par,var)
if not (par::numeric and var::numeric) then
return 'procname'(args);
end if;
if paramloc <> par then
paramloc := par;
WV(parameters=[U0=paramloc]);
end if;
WV(var);
end proc:
end module:

For fun, this time I'll use the parameter calling sequence to flip the
axes, instead of plots:transform. That's just because I want t displayed
on the first axis. But for the performance gain, what matters is that it
is U0 which gets values from the first axis plottingrange.
> 
str := time[real]():
plot3d([y,x,VofU0module(x,y)], x=U0lo..U0hi, y=tlo..thi,
grid=[n,m], labels=["t","U0","V(t)"]);
(time[real]()str)*'seconds';

And, naturally, I could also use the parametric form to get a fast plot
with the axes roles switched.
> 
str := time[real]():
plot3d([x,y,VofU0module(x,y)], x=U0lo..U0hi, y=tlo..thi,
grid=[n,m], labels=["U0","t","V(t)"]);
(time[real]()str)*'seconds';

