A few issues and questions.

!. Most of us prefer not to have to type your code. So bring the code as text, please.
2. Don't you intend the end value of the loop parameter to be p-1 rather than 1-p?
3. What is x=x intended to do?
4. What does disp do?
5. Your loop will run to the end and x will finally have the value loop_end+1 (i.e. p if I'm right in point 2).
6. Presumably you want something like this:
   if (a^x-b) mod p<>0 then return x elif (a^x-b) mod p = 0 then return disp(x) end if
7. What role does z play?
8. Why is x declared global?

@Preben Alsholm If we keep the 8 boundary values as you had them, although simplifying them to

bcs11 := f(0) = S, D(f)(0) = 1, D(f)(blt) = 0, D(F)(blt) = 0, F(blt) = f(blt), theta(0) = 1, theta(blt) = 0, theta1(blt) = 0;

and replace Eq2 by
Eq2d:=diff(Eq2,eta);
#then with gamma replaced by 1:
res:=dsolve(eval({Eq1,Eq2d,Eq3,Eq4,bcs11},{A=L[1],gamma=1}),numeric,initmesh=256);
plots:-odeplot(res,[[eta,f(eta)],[eta,F(eta)],[eta,theta(eta)],[eta,theta1(eta)]]);
#Now the question is: How well is the original Eq2 satisfied?
#It turns out not to be too bad at all:

plots:-odeplot(res,[eta,eval(lhs(Eq2),A=L[1])]);

That should not be a surprise for the fpllowing reason:
Because of Eq2d we have that the solution to the new problem Eq1,Eq2d,Eq3,Eq4,bcs11 satisfies lhs(Eq2) = constant. But since
subs(eta=blt,bcs11,convert(Eq2,D));
returns 0=0 we see that that constant is 0.
######
I may add that this same result is also obtained by using rifsimp with the dependent variables specified as in this command:

T:=DEtools[rifsimp](eval([Eq1,Eq2,Eq3,Eq4],{gamma=gamma1}),[theta,theta1,F,f]):
#This time there is a constraint equation just as we had Eq2 above:
indices(T);
T[Pivots];
T[Case];
T[Constraint];
indets(T[Solved],specfunc(anything,diff));

I think the crucial difference is shown immediately at the evaluation of the calls to sqrt:

restart;
S:=sqrt(4-sqrt(7));
S1:=sqrt(4-sqrt(6));

##sqrt is a module that exports ModuleApply, so sqrt can be used as a procedure sqrt:
type(sqrt,`module`);
eval(sqrt);
showstat(sqrt);

What is the actual function?

That was a retorical question.

BesselJ(0,k) is an actual function. For a definition go to
?BesselJ

To get an infinite series convergent for all k try

convert(BesselJ(0,k), Sum);
##It looks nicer this way:

subs(_k1=n,convert(%,factorial));

#######
I suppose that what you actually were working with was piecewise(abs(x)<1,1/sqrt(1-x^2),0) and not 1/sqrt(1-x^2):

inttrans[fourier](piecewise(abs(x)<1,1/sqrt(1-x^2),0),x,k);

I'm just using the pasted code lines above. However, you should be solving for x and y. If you don't give solve a second argument it assumes that you want to solve for all names (including phi and theta).

To get started try:
rrestart;
eq1 := x*(1.6*(1-(1/100)*x)-phi*y);
eq2 := (x/(15+x)-0.3e-1*x-.4)*y+.6+theta;
solve({eq1,eq2},{x,y},explicit);
de1 := diff(x(t), t) = x(t)*(1.6*(1-(1/100)*x(t))-phi*y(t));
de2 := diff(y(t), t) = (x(t)/(15+x(t))-0.3e-1*x(t)-.4)*y(t)+.6+theta;

res:=dsolve({de1,de2,x(0) = 990.0, y(0) = 10},numeric,parameters=[phi,theta]);
res(parameters=[phi=0.5,theta=5]);
plots:-odeplot(res,[[t,x(t)],[t,100*y(t)]],0..10);

## About your worksheet:  Do not define phi or theta as arrays, keep them as names!

There are a few problems.

1. Did you intend gamma in Eq4 to be Euler's constant? If not use another name like gamma1 and assign a value to gamma1.

2. The procedure `dsolve/numeric/bvp/convertsys` cannot isolate the highest order derivatives uniquely, which is why you get your error statement. Try
solve({Eq1,Eq2,Eq3,Eq4},{diff(f(eta),eta$3),diff(theta(eta),eta,eta),diff(F(eta),eta),diff(theta1(eta),eta)}):
indets(%,RootOf);
You will see that a solution is obtained, but involves an unsolved quadratic in _Z. This is simply because the highest order derivative for F is F' and that appears squared in Eq2.
To find the highest orders you can look at the output from
indets({Eq1,Eq2,Eq3,Eq4},specfunc(anything,diff));

3. Your boundary conditions also have problems. There are 8, but the total order of the system is 7:
[F: 1, f: 3, theta: 2, theta1: 1] thus total order =1+3+2+1 = 7.
Furthermore, you cannot have boundary conditions involving derivatives of the highest order (or above). You have D(F) in your bcs1.

To get over the problem mentioned in point 2 you could replace the RootOf with one of the two roots, but which?
############
## Trying that idea with gamma1 = 1 and with corrected boundary conditions:
bcs1a:=f(0) = .5, D(f)(0) = 1, D(f)(5) = 0, F(5) = f(5)+5*D(f)(5), theta(0) = 1, theta(5) = 0, theta1(5) = 0;
#Taking the first of the two roots when solving Eq2 for diff(F(eta),eta):
Eq2a:=diff(F(eta),eta)=solve(Eq2,diff(F(eta),eta))[1];
res:=dsolve(eval({Eq1,Eq2a,Eq3,Eq4},{A=L[1],gamma=1}) union {bcs1a},numeric);
odeplot(res,[eta,f(eta)]);
#That actually put out a result!
############

Another idea is to try DEtools[rifsimp]. This procedure doesn't only simplify algebraically but also differentiates. When successful the returned system should be equivalent to the given one.
T:=DEtools[rifsimp]([Eq1,Eq2,Eq3,Eq4]):
The output is large, so I have suppressed printing it.
T will be a table, try:
indices(T);
You will see [Solved], [Pivots], [Case].
These entries are fairly short:
T[Pivots]; #Reminding you of gaussian elimination
T[Case]; #Assumptions for isolating the highest derivatives
The most interesting is T[Solved], which is large. It contains the new system. You can try
indets(T[Solved],specfunc(anything,diff));
Compare with the similar command for your original system.
Find the left hand sides:
lhs~(T[Solved]);
has(rhs~(T[Solved]),%);
The outputs from these two commands show that the highest derivatives have been isolated.

I tried:
restart;
C:=convert(VectorCalculus:-PositionVector([u,v,w], spherical),list);

## and got  C := [u*sin(v)*cos(w), u*sin(v)*sin(w), u*cos(v)]

## Then I looked up:

?coords
##and found in the text the definition of sperical coordinates:
                        x = u cos(v) sin(w)
                        y = u sin(v) sin(w)
                            z = u cos(w)
## Now animate the drawing of a unit sphere in two ways:
C1:=subs(u=1,C);
plots:-animate(plot3d,[C1,v=0..Pi,w=0..pp,scaling=constrained],pp=0..2*Pi);
plots:-animate(plot3d,[C1,v=0..pp,w=0..2*Pi,scaling=constrained],pp=0..Pi);
##Now look up
?spherical coordinates
## Here I find the example
plot3d(1,t=0..2*Pi,p=0..Pi,coords=spherical,scaling=constrained);
## which I animated in two versions:
plots:-animate(plot3d,[1,t=0..2*Pi,p=0..pp,coords=spherical,scaling=constrained],pp=0..Pi);
plots:-animate(plot3d,[1,t=0..pp,p=0..Pi,coords=spherical,scaling=constrained],pp=0..2*Pi);

Yes, I am also confused.

You certainly have a problem starting with tm:

tm := t->t  mod 86400;
tm(6.1); #tm only accepts integers

Thus all the functions after that building on tm suffers from the same problem.
When you define f__0 (the strange tau is an artifact of the detestable 2D input):
f__0 := z-> int(1/`&tau;__0`(x), x = 0 .. z);*
you have a genuine problem since the definition of `&tau;__0` relies on tm.
Try
`&tau;__0` (6.1);
Thus the integral makes no sense.

If only positive values of x or z are considered you could consider changing the definition of tm to

tm:=t->86400*frac(t/86400);
and the definition of f__0 to
f__0 := z-> evalf(Int(1/`&tau;__0`(x), x = 0 .. z));

This latter assignment makes use of numerical integration (evalf(Int(...)) instead of int). The following actually turns out a result
f__0(20);
namely, 0.0005577834770.

The way you have defined eqns, ics, var, data, they are all sequences, not sets nor lists.
Keeping eqns, ics, var, and data as you have defined them you can just do as follows:

res := dsolve({eqns, ics}, {var}, method = laplace);

where I have used method=laplace which comes with an answer very fast. It involves sums over the roots of a quartic polynomial, which is to be expected as your system consists of 4 first order equations.

To continue with inserting parameters, do

eval(res,{data}); #Notice {} again
allvalues(%);
evalc(%);
RES:=evalf(%);
Z:=subs(RES,[var]); #Notice list brackets []
plot(Z,t=0..5);

##Now inserting parameters before solution:
resC := dsolve(eval({eqns, ics},{data}), {var}, method = laplace);
RESC:=evalf(evalc(allvalues(resC)));
ZC:=subs(RESC,[var]); ##Notice list brackets []
plot(ZC,t=0..5);
##One should expect the two list Z and ZC to agree except for errors caused by using floats.
Z-ZC;

List brackets were used in producing Z and ZC in order to be certain of the order.

Incidentally, no use is being made of your assumptions. You can safely remove them.
##########
A little linear algebra often comes in handy:
LinearAlgebra:-GenerateMatrix(rhs~([eqns]),[var]);
A:=%[1];
pol:=LinearAlgebra:-CharacteristicPolynomial(A,lambda);
polC:=eval(pol,{data});
solve(polC,lambda);
evalc([%]);
evalf(%);

ODE:= diff(y(x),x$2) - 1.0325*diff(y(x),x) + 1.36*y(x) = sin(2*x):
ODE3:=diff(ODE,x);
bcs:={y(0)=0, y(1)=1};
bcs0:=subs(x=0,bcs,convert(ODE,D));
S:=dsolve({ODE3,bcs0} union bcs,numeric);
plots:-odeplot(S,[x,diff(y(x),x,x)],0..1);

In the help page for dfieldplot it doesn't mention initial conditions as acceptable. On the contrary it says,

"For plotting solution curves, see DEplot or phaseportrait."

Now your system is not autonomous because of the piecewise defined de1.
However, both of these are autonomous:

de11:=simplify(de1) assuming t<8;
de12:=simplify(de1) assuming t>8;
dfieldplot([de11,de2], [A(t),G(t)],t=0..1,A=0..900,G=0..200);
dfieldplot([de12,de2], [A(t),G(t)],t=0..1,A=0..900,G=0..200);

Notice that the t-range t=0..1 might as well be t=90..6789 (or whatever). It is not used, but has to be there.

If I were you I would only use DEplot. You can do everything with that.

To find nullclines you don't solve eq1 and eq2 together. You do one at a time:

solve(eq1,{u(t),v(t)}); #Finding the vertical nullclines
solve(eq2,{u(t),v(t)}); #Finding the horizontal
V:=plot([[0,v,v=0..1.2],[2-2*v,v,v=0..1.2]],color=blue,thickness=3,view=[0..2,0..1.2]);
H:=plot([0,[1,v,v=0..1.2]],color=green,thickness=3,view=[0..2,0..1.2]);
#Suppose the phaseportrait has been assigned to the variable ph, then do:
plots:-display(ph,H,V);

Just do:
H,H1:=op(subs(S1,[y(x),diff(y(x),x)]));

plot([H(t),H1(t)], t = 0..1, thickness = 4);

#For plotting purposes you don't need to define H and H1:
plots:-odeplot(S1,[[x,y(x)],[x,diff(y(x),x)]],0..1);

As in the problem in your question a week or two ago this problem has symmetry about x=0.
Using a similar approach we set
sol:={phi(x)=A1*cosh(r1*x)+B1*cosh(r2*x)+C1*cosh(r3*x),
     psi(x)=A2*cosh(r1*x)+B2*cosh(r2*x)+C2*cosh(r3*x)};

Here r1,r2,r3 are the 3 positive roots of the characteristic polynomial for the matrix A obtained by rewriting your system as a first order system in the traditional way. We must expect A2,B2,C2 to depend on A1, B1, C1.

The characteristic polynomial for A is

pol:=lambda^6+(a11*a77+a22*a66-a44^2)*lambda^4/(a11*a66)+(a22*a77+a33*a66-2*a44*a55)*lambda^2/(a66*a11)+(a33*a77-a55^2)/(a66*a11);

With your constants and Digits:=20 you get the following positive roots:

r1, r2, r3 := .57466917365728616069, .69502623409256130812, 1.2578404958328309236;

Now insert sol into eq1 using odetest:
odetest(sol,eq1);
res:=collect(%,cosh);
Eqs:=[seq(coeff(res,cosh(r*x)),r=[r1,r2,r3])];
solve(Eqs,{A2,B2,C2});
sol2:=eval(sol,%);
eval(eq2,sol2); #Just a test, but it ought to be exactly zero, however rounding ...
#Now to determine A1,B1,C1 we use the boundary conditions on the right:
bcs1:={phi(a)=sigma1,D(phi)(a)=0,psi(a)=sigma2};
subs(sol2,phi(x));
Eqsb1:={eval(diff(%,x),x=a)=0}; #The derivative condition for phi
Eqsb2:=subs(x=a,bcs1,sol2); #The two other conditions
fsolve(Eqsb1 union Eqsb2,{A1,B1,C1});
#Answer: {A1 = 3.8652147206626085692*10^7, B1 = -1.9405216687792465984*10^7, C1 = -1.0920477631668062160*10^6}
#However, continue with
sol3:=eval(sol2,%);
which is the solution.
plot(subs(sol3,[phi(x),psi(x)]),x=-a..a);

Try this version

f := proc(kValue, KAValue) global ODE_Solution, k, K, C, t;
  ODE_Solution('parameters' = [k = kValue, K[A]= KAValue]);
  try
    return eval(C[A](t), ODE_Solution(input));
  catch:
    return 100;
  end try;
end proc;

1. The construction a[b]d is not an acceptable name. You use such a construction as one of your formal parameters.
2. K[A] refers to the entry of the table K with index A. So the global declaration should involve K, not K[A].