@Naeem Ullah You have Euler's constant gamma ( evalf(gamma) = .5772156649) in some places and gamma1 in other, where gamma1 is later defined to be 1. Is this intentional or should both be gamma1?

It is rather strange to me to see abserr=1e20 considering that the default is 1e-6.

Clearly you don't use or need the deprecated linalg package. If you need LinearAlgebra somewhere later, use LinearAlgebra.

Do you really want to set Digits:=24?

@nm As one user who has answered (or tried to answer) an awful lot of questions about boundary value problems for third or fourth order odes I must agree with Carl; there have been very many, but most have had to do with convergence problems, not with the syntax for higher derivatives in the bcs.
Looking under my own name in Users, All answers, first page, the first question about a bcs problem is this one in which the second derivative happens to appear in the bcs:
http://www.mapleprimes.com/questions/218225-Error-Showing-As-#answer232183

@artfin I wasn't referring to your use of t->...., but to the first argument to Int, which wasn't a procedure, but was a procedure call, i.e. sin(x) versus sin.

The t from the error doesn't come from the obvious t in the outer construction of yours:
comp_1:=t->evalf(Int(x->r1_lab_deriv(x)[1],0..t,epsilon=1e-6,method=_d01akc)):

as will be seen if you change that t to e.g. s (in both places). The error still refers to t.

The error is due to the t in r1_lab_deriv and inside that it is due to the term r1_rot_deriv(t).
r1_rot_deriv(t) you have defined as
r1_rot_deriv:=x-> <fdiff(r1_rot(t)[1],t=x),fdiff(r1_rot(t)[2],t=x),fdiff(r1_rot(t)[3],t=x)>: using fdiff.
Furthermore
r1_rot(t) is defined as
<-r1(t)*cos(q(t)/2),0,-r1(t)*sin(q(t)/2)>:
where r1 and q are procedures returned by dsolve/numeric (result in solution_procedure).
The error comes from the fact that procedures returned by dsolve/numeric with output listprocedure return unevaluated when receiving just a name like t.
Since it follows from eqs that diff(r1(t), t) = 0.5443205700e-3*p1(t) you could define
r1_deriv by 0.5443205700e-3*p1 and similarly q_deriv by  using the odes.
The less convoluted the end procedure is the better. You have sequences of procedures calling each other.
##
But it may be better (and easier) to include needed procedures in the original call to dsolve/numeric.
I shall confine myself to a simple example:
##
restart;
ode:=diff(x(t),t)=sin(x(t))+t;
res:=dsolve({ode,x(0)=0},numeric,output=listprocedure);
X:=subs(res,x(t));
fdiff(X,[1],[1.2345]);
eval[recurse](rhs(ode),{t=1.2345,x(t)=X(t)});
## Now include the derivative as xp instead:
res2:=dsolve({xp(t)=rhs(ode),lhs(ode)=xp(t),x(0)=0},numeric,output=listprocedure);
X,XP:=op(subs(res2,[x(t),xp(t)]));
XP(1.2345);

@Naeem Ullah There is no way we can help unless we get the entire code, either as text (not images) or as an uoloaded worksheet. Use the fat green arrow in the second line to the right in the MaplePrimes editor.
We need to be able to work with the code ourselves without having to type in your code first.

@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.