What has this got to do with dsolve?

@Christian Wolinski 

In Maple 8 the command
 type('u[1](x)', specfunc(anything, u));

returns false.

In Maple 12 (and Maple 2016) it returns true.

For your other command
u:=proc(x) x+1 end; v:=table([]); type('u[1](x)', specfunc(anything, u)), type('v[1](x)', specfunc(anything, v));

Maple 8 returns false,false. Maple 12 and Maple 2016 return true,true.

You write:
"I would like to differentiate the (Sst/Vst)/(Spt/Vpt) by dsm first then integrate it with dsm ranges from 0 to dmax to get my final answer...".

That seems to me to correspond to f(x) -f(0) = int( diff(f(t),t), t=0..x);

Why not just do:
u:=Sst*Vpt/(Vst*Spt);
eval(u,dsm=dmax)-eval(u,dsm=0);

Helmut Weber has provided C++ source code for a modification of Talbot's algorithm due to L.N. Trefethen, J.A.C. Weideman, and T. Schmelzer in Talbot quadratures and rational approximations, BIT 46, 653-670  (2006).

http://www.cs.hs-rm.de/~weber/lapinv/talbot/talbot.htm

This code is easily rewritten in Maple and is very fast. The behavior already shown above is confirmed: The inverse laplace transform f(t) of (5*s^(3/10)+s)/(s^2+5*s^(13/10)+5*Pi^2) dies off pretty rapidly after t=3. It seems to be creeping up to 0 from below as t->infinity.

Trying the method sketched above with monomials in t (t^n, n=0..10), but this time only on the interval t=0..3 and sampling at mapped chebyshev abscissae I got
f(t) = 0.1679667698e-1*t^10-.2316408948*t^9+1.173311393*t^8-1.958556676*t^7-4.80665368*t^6+
28.70316704*t^5-56.34317461*t^4+52.98096151*t^3-20.63139377*t^2+1.
which plots like this:

I compared that with a numerical inversion of the laplace transform found earlier to be
(5*s^(3/10)+s)/(s^2+5*s^(13/10)+5*Pi^2)

I used a contour integral for inversion using a contour apparently due to Talbot. I found that in:
http://math.arizona.edu/~brio/NLAPUNR.pdf

See also
http://math.arizona.edu/~brio/WEEKS_METHOD_PAGE/pkanoTalbotsMethod.html

I used sigma=0, mu=1, nu=3 in the contour given by s = sigma+mu*theta*(I*nu+cot(theta)).
I computed f(t) at t-values from 0.01 and up until I ran into a problem at t = 2.80.
The resulting graph was almost indistinguishable from the one presented here.

@farahnaz The names _C1, _C2, etc. are used by pdsolve (and dsolve) as names of arbitrary constants.

In your case this would imply that you should determine the constants so that sol3[3] satisfies the requirements for the behavior in rho at 0 and infinity. That may be very difficult or practically impossible.

By using the laplace transform or by separating variables it appears that the solution can be written
T(z,t) = sin(Pi*z)*f(t)

where f(t) satisfies the following integral equation:
eq:=simplify(eval(PDE,{T=((z,t)->sin(Pi*z)*f(t))})/sin(Pi*z));
    
         

The initial values must be:
ics:={f(0)=1,D(f)(0)=0}; #Based on the initial values for T(z,t)

I tried with my own unsophisticated solver to find an approximate solution to eq with ics first in terms of monomials in t and on the interval 0..10. Since that produced ever increasing waves when animating T(z,t) I turned to using cos(n*Pi*t), n=0..N instead of the monomials.
I shall return if I get anything I believe could be close to an acceptable approximate solution.

@farahnaz When changing from dsolve to pdsolve I get what I (to great length) described.
As an example sol3[1] is found to be

H(rho, z) = (1/2)*Pi*(_C1*BesselJ(1, rho)+_C2*StruveH(1, rho))/_C2

which (as I also wrote) is independent of z.

I just edited my answer to include a correction: sol3[1] doesn't satisfy the requirement: I had them turned around.

Since sol3[2] doesn't either, you are stuck with sol3[3] which is the complicated one.

@tomleslie I just tried your code in Maple 12 after having replaced all occurrences of numelems with nops.
It worked nicely!
So surely it will work in Maple 13 too.

@shakuntala You don't update from Maple 13 to a newer version: You simply buy another version, so buy Maple 2016 since money is not a problem. Go to the Maplesoft website for purchase.

@torabi The command is what you used and what I used in the form
invlaplace(rhs(sol3),s,t);

It was not successful, however. invlaplace couldn't do the inversion, but just used linearity to get

Why is the title "Shooting Method" ?
There is no shooting going on. `dsolve/numeric/bvp`doesn't use a shooting method.

@Bendesarts Using your condition1 I tried

indets(condition1,`^` );

There are no squared sine functions among the members of the output:

@Christopher2222 With few exceptions (e.g. Edgardo S. Cheb-Terrab and Robert J. Lopez), Maple employees rarely participate on MaplePrimes, but leave the discussion to the members.
I'm sure there are very good reasons for that, among those

(1) too busy working on matters directly related to Maple software;
(2) allowing a free discussion among members without interfering with an "authoritative" or "official" verdict on the issue.

Both of those two are perfectly legitimate and honorable reasons.

But the spam issue is different as it has nothing to do with the behavior of Maple software, so some reaction from the responsible part would be appreciated.

@9009134 Why not keep the RootOf. Notice that it is indexed with 1,2,and 3.
Try e.g. setting all parameters to 1:

nms:=indets(A,name);
eval([Sol],nms=~1);
allvalues(%);
evalf(%); #If you want floats you can skip the previous step.

I didn't even think of trying to integrate u1 symbolically before, but doing that improves the result.
For convenience I bring all of the code.
I have made two extra changes. In series I have used the order 5 instead of 3.
Also it seems to help to use the procedural form for input to evalf/Int:

restart;
Digits:=15:
tD := 20;
u:=4/Pi^2*( 1 - exp(-x^2*tD) ) / ( x^3*( BesselJ(0,x)^2 + BesselY(0,x)^2 ) );
series(u,x,5); #Replaced 3 by 5
u1:=convert(%,polynom);
U1:=int(u1,x=0..1/2);
#evalf(Int(u-u1,x=0..1/2,epsilon=1e-6));
evalf(Int(unapply(u-u1,x), 0..1/2, epsilon=1e-12));
evalf(U1);
res1:=%%+%;
res2:=evalf(Int(u,x=1/2..infinity));
res1+res2;
simplify(%); # Result12.3201290912317