@9009134 I really don't know if there is a solution or not.
If you have any idea how a possible solution for F, K, Omega, and theta might look then I would like to know.
That could be used in dsolve/numeric/bvp as an approximate solution.
The input would be of the form 
approxsoln=[F(eta)=... , K(eta) = ... , Omega(eta) = ... , theta(eta) = ...] 

where of course the dots should be replaced with an algebraic expression in eta.

I have only been able to come up with results that mostly look like F(eta) = K(eta) = Omega(eta) = 0 for all eta=0..1 (actually just very small, but that worries me) and theta(eta) = 1-eta.

I can show you a plot of  [10^3*F,10^6*K,10^8*Omega,theta] where it is important to notice the scale factors used in order to fit all 4 graphs into one. I should remark too that I obtained this from dsolve/numeric/bvp with input from a private program (written by myself) used as an approximate solution, but I only achieved convergence from dsolve with abserr=5e-3.
If these graphs look promising, let me know.

@Federiko I tried replacing f by a linear combination of s^n, n=0..10 (coeffs c(n) ) and sampled the result at Chebyshev zeros, after which the integration was performed and solved for the coefficients c(n).

What I got was:
res := f(s) = (115.4500366-40.68349123*I)*s^10+(-455.2724280+183.4086956*I)*s^9+(546.1618683-333.2102196*I)*s^8+(138.6309983+293.1462345*I)*s^7+(-882.5976097-103.8372237*I)*s^6+(711.7711146-6.166277135*I)*s^5+(-96.85736837-6.915592829*I)*s^4+(-98.03396901+19.04906798*I)*s^3+(-.8642921380-0.56852810e-1*I)*s^2+(21.25808967-5.232446726*I)*s-0.9123647e-4-0.63278e-5*I

which is not terribly good (and not terribly bad) considering:
test:=eval((lhs-rhs)(eq),f=unapply(rhs(res),s)); #where eq is your integral equation.
numapprox:-infnorm(abs(test),s=0..1); # result 0.001435744964
## Plotting is faster:
plot(abs(test),s=0..1);

Maple doesn't do integral equations numerically. For which interval of s-values do you want to find f(s)?

I had a look at your system, which I shall call SYS (including the boundary conditions.
Obviously you are looking for a solution (F,K,Omega,theta) satisfying K(eta)<>0, Omega(eta)<>0 for all eta in the open interval (0,1).
Whether such a solution exists is not clear to me. It may not.
It may be noticed that theta only appears in its "own" differential equation, i.e. you can split the system into 3 odes and boundary conditions for F,K,Omega having no theta in them at all.
All the boundary conditions for F,K,Omega have zero values.

restart;
SYS:={diff(theta(eta), eta, eta)-3*Omega(eta)*(F(eta)*(diff(theta(eta), eta))-theta(eta)*(diff(F(eta), eta)))/(2*K(eta))+((diff(K(eta), eta))/K(eta)-(diff(Omega(eta), eta))/Omega(eta))*(diff(theta(eta), eta)) = 0, diff(F(eta), eta, eta, eta)+Omega(eta)*(3*F(eta)*(diff(F(eta), eta, eta))-(diff(F(eta), eta))^2)/(2*K(eta))+((diff(K(eta), eta))/K(eta)-(diff(Omega(eta), eta))/Omega(eta))*(diff(F(eta), eta, eta))+Omega(eta)/K(eta) = 0, diff(K(eta), eta, eta)+Omega(eta)*(1.5*F(eta)*(diff(K(eta), eta))-K(eta)*(diff(F(eta), eta)))/K(eta)+((diff(K(eta), eta))/K(eta)-(diff(Omega(eta), eta))/Omega(eta))*(diff(K(eta), eta))+(diff(F(eta), eta, eta))^2-Omega(eta)^2 = 0, diff(Omega(eta), eta, eta)+Omega(eta)*(3*F(eta)*(diff(Omega(eta), eta))+Omega(eta)*(diff(F(eta), eta)))/(2*K(eta))+((diff(K(eta), eta))/K(eta)-(diff(Omega(eta), eta))/Omega(eta))*(diff(Omega(eta), eta))+Omega(eta)*(diff(F(eta), eta, eta))^2/K(eta)-Omega(eta)^3/K(eta) = 0, F(0) = 0, K(0) = 0, Omega(0) = 0., theta(0) = 1, theta(1) = 0, (D(F))(0) = 0, (D(K))(1) = 0, (D(Omega))(1) = 0, ((D@@2)(F))(1) = 0};
### Now splitting as described above:
sys_theta,sys_FKO:=selectremove(has,SYS,theta);
## Clearly, if you can solve the original problem SYS then you can also solve sys_FKO and if you cannot solve sys_FKO there is no chance for SYS.
##
The bad news: So far what I have been able to come up with is not trustworthy.

@acer In my experience a lot of Danes have had this problem. It could be because Danes are rather stupid (unlikely, I'm one), or that special characters æ,ø,å cause problems. But there are quite a few other languages using special characters.
So what is the reason?

I learned years ago (from bad experiences with my students' worksheets) never to save a worksheet with output. Always delete the output before saving.
I still don't save output ever. These days I never have a problem although I use a Danish keyboard.

@Christian Wolinski No, same result in all versions I believe:

type(f[1](x),indexed); # false
But of course you can do:
type(op(0,f[1](x)),indexed); # true

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.