For the fun of it (seriously, yes) you could try deactivating the error check by setting abserr ridiculously high.
It won't hurt to try. Here I'm using abserr=10^10, which basically means that that all components of the corresponding first order system are allowed to stray 10^10 away from the "true" solution!!!
You won't be allowed to use abserr = infinity, but 10^10 is just a substitute.

Digits:=15:
res3 := dsolve(sys2 union bcs3, numeric, maxmesh = 2048, abserr = 10^10,initmesh=256);
plots:-odeplot(res3, [[y,g1(y)],[y,100*g2(y)]], 0 .. 1);


Whether this has anything whatsoever to do with a solution I don't know.
As a warning you may try plotting the derivatives of g2. Here the third (corrected from 2nd) derivative:
plots:-odeplot(res300, [y,diff(g2(y),y,y)], 0 .. 1);

@9009134 I forgot to give the plotting command that produces the graph shown:

plots:-odeplot(solA,[seq([eta,fu(eta)],fu=[10^3*F,10^6*K,10^8*Omega,theta])],0..1,legend=[10^3*F,10^6*K,10^8*Omega,theta],color=[red,blue,green,black]);

As I said I got inspiration to the approximate solution from results obtained with my private program. These results were in matrix form. They were used as an approxsoln in dsolve/numeric/bvp and the graphs in my comment (the first set) were produced. Then second time around I simply looked at those graphs and took out of my head the approxsoln, which basically was
[F(eta)=eta^2,K(eta)=tanh(eta),Omega(eta)=tanh(eta),theta(eta)=1-eta]
but subsequently modified with scaling factors to make the graphs of those look like the graphs in my comment.
Thus I wasn't doing anything fancy this second time.

It is disturbing that the second set of graphs (the graphs in the "answer") although having the same form are clearly different from the set of graphs in my "comment", e.g. does the graph of 10^6*K level off at about 2.5 in the "answer", but at about 1 in the "comment". Nonuniqueness is possible with your boundary conditions, but I wouldn't bet much on the results.
Please see,however, the addition to my answer where I replace the requirement (D@@2)(F)(1)=0 by F(1) = 1.
##
In the help page for dsolve/numeric/bvp it also clearly states that the methods used are not well suited for problems with singularities.

@Naeem Ullah I don't see a link to the worksheet in your latest reply.

@Naeem Ullah But you need to give that gamma2 a value. What would it be? The value will obviously influence the results.

By a greater value for abserr I simply meant a value greater than 1e-6, as e.g.  1e-3 .
You may even go back to "inactivating" the error check by setting it ridiculously high like you had originally (10^20) just to see what you get.

It is worth mentioning that the help page for dsolve/numeric/bvp says this about the methods used:
... the midpoint submethods are capable of handling harmless end-point singularities ...

and continues with

The available methods are fairly general, and should work on a variety of BVPs. They are not well suited to solve BVPs that are stiff or have solutions with singularities in their higher order derivatives.

As for contact through email I prefer to discuss the issues in this forum.


 

@Naeem Ullah You said that gamma is the curvature parameter and gamma1 is just a constant set at 1.

But gamma in your new version of the worksheet is still just Euler's constant, and is an initially known constant in Maple.
Thus it couldn't be a curvature parameter as it appears now. It will be treated as a constant with value approximately  .5772156649 as I said earlier.

You never assign a value to gamma, and in fact you won't even be allowed to:
gamma:=1.2345; # results in the error:

Error, attempting to assign to `gamma` which is protected.  Try declaring `local gamma`; see ?protect for details.

As the error message says: You can start your worksheet like this:
restart;
local gamma;
##Then the rest. Now you can assign to gamma.
I find it simpler just to use another name instead of gamma, e.g. gamma2.

About abserr: I didn't mean to say that you have to use the default abserr, which is 1e-6. You may have to raise it.
I was just pointing out that abserr=10^20 seems ridiculous as it would mean that all that is guaranteed about the solution is that it hopefully would not be more than 10^20 away from the true solution!
It does appear, however, that dsolve/numeric/bvp runs through some initial work which may bring forth a solution that is actually not bad, but you will have no guarantee that it is any good.
We may conclude that using abserr=10^20 or some huge number like that basically makes dsolve/numeric/bvp skip error checking.
Skipping error checking isn't all that bad an idea if you are desperately trying to get some idea of how the solution looks.

@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