The operator Mf of multiplication by the function f can be defined like this:

Mf:=g->f*g;
#Applying Mf to a function h (in the mathematical sense, not Maple)
(Mf(h))(x); #The parenthesis around Mf(h) i.e. (Mf(h)) are superfluous, but are added for clarity of meaning.
((D-Mf)@@2)(g)(x);
expand(%);

#An operator M which maps functions onto multiplication operators:

M:=f->(g->f*g);
#Thus M(f) is the operator which acts like this:
M(h)(g)(x);
((D-M(f))@@2)(g)(x);
expand(%);

########Multiplying by a constant function:
M(89)(g)(x);
a:=x->k; #The function a defined by a(x)=k for all x
M(a)(g)(x);

############## Followup question
#To accomplish your intended result you can use M above:

restart;
M:=f->(g->f*g);
id:=(x,m)->x;
M(id)(g)(y,n); #test
(M(id)@D[1])(g)(x,m); #test
((M(id)@D[1])@@2)(F)(x,m);

#Alternatively use id but not M:
Mx:=g->id*g;
((Mx@D[1])@@2)(F)(x,m);

You can use approxsoln (after correcting the B-problem pointed out by Carl):

AP:=NULL:
for k to 6 do
R := dsolve(eval({Eq1, Eq2, bcs1, bcs2}, B = L[k]), [f(eta), theta(eta)], numeric, output = listprocedure,maxmesh=512,AP);
AP:=approxsoln=R;
X1 || k := rhs(R[3]); X2 || k := rhs(R[4]); Y1 || k := rhs(R[5]); Y2 || k := -rhs(R[6])
end do;

The print command at the end then gives me the values:

[-1.41421356296431, -1.43879570677334, -1.46322186668648, -1.48747814527757, -1.51155337283017, -1.55912751987572]

Incidentally, there is no reason for the 'print'. Just writing
[(X2 || (1 .. 6))(0)];

will do.

I'm answering your question about tangency points.

Since fsolve is rather reluctant to return a result if it not perfect, or if it doesn't appear so to fsolve, I once wrote a procedure to pick the best solution found by fsolve in a failed attempt.
It is given below and called BestPoint.
It requires as input either a procedure or a list of procedures with option remember.
It finds the input argument(s) corresponding to the least value of an error function applied to the output of the procedure(s).
To test it just copy and paste the whole thing and execute it in the order given.
I'm applying it to your tangency problem.
The procedure on which BestPoint will be used is called q below.

restart;
BestPoint:=proc(q::{procedure,list(procedure)},{errorfunction::procedure:='default',procedurelist::truefalse:=false,
number::posint:=`if`(not procedurelist and q::list(procedure),nops(q),1)}) local Q,f,x,i,T,T1,E,R,S,inf;
 inf:=proc(eq) lhs(eq)=evalindets~(rhs(eq),Not(numeric),x->infinity) end proc;
 if errorfunction ='default' then f:=unapply(add(x[i]^2,i=1..number),seq(x[i],i=1..number)) else f:=errorfunction end if;
 if not procedurelist then
    if not type(q,list) then Q:=[q] else Q:=q end if;
    T:=op~(4,eval~(Q)); #Remember tables
    if T=[] then error "No remember tables" end if;
    T:=table~(map(inf~,op~(eval~(T)))); #replacing nonnumeric entries with infinity
    E:=[entries]~(T,nolist);
    R:=ListTools:-Enumerate(f~(op(E)));
    T1:=T[1]    
 else
    T:=op(4,eval(q));
    if T=NULL then error "No remember table" end if;
    T:=table(map(inf,op(eval(T)))); #replacing nonnumeric entries with infinity
    E:=[entries](T,nolist);
    if number<>nops(E[1]) then error "Wrong or missing argument: number = n" end if;
    R:=ListTools:-Enumerate((f@op)~(E));
    T1:=T
 end if;
 S:=sort(R,(x,y)->abs(x[2])<abs(y[2]));
 print(ErrorFunction=eval(f),Error=S[1,2]);
 indices(T1)[op([1,1],S)]
end proc:

#Your tangency problem.
Output from q is a list of two real numbers. These numbers we want to make as small as possible and will use fsolve on q in order to do that. This will fail in many cases, but then we use BestPoint.
On my run of the H-loop at the bottom (with Digits=10) fsolve succeeded 5 times out of 14 in Maple 2015.

q := proc (x,y,H) option remember;
local f0,f1, f2, k;
  f0:= 0.585=abs(k*(0.585+H*cos(20*Pi/180))+0.585-H*sin(20*Pi/180))/sqrt(k^2+1);
  f1:=fsolve(f0,k)*x-y;
  f2:= (x-0.585-H*cos(20*Pi/180))^2+(y+0.585-H*sin(20*Pi/180))^2-0.585^2;
  evalf([f1,f2])
end proc;
q1:=proc(x,y) q(x,y,H)[1] end proc;
q2:=proc(x,y) q(x,y,H)[2] end proc;
#fsolve cannot take q as input; that is the basic reason for q1,q2.

H:=-0.18: sol := fsolve([q1,q2]);
res:=table();
res[1]:=BestPoint(q,procedurelist,number=2);
printlevel:=2:
i:=1:
for H from -0.18 by 0.03 to 0.18 do
  i:=i+1;
  sol:=fsolve([q1,q2],res[i-1][1..2]);
  if not type(sol,list(numeric)) then res[i]:=BestPoint(q,procedurelist,number=2) else res[i]:=[op(sol),H] end if;
  forget(q)
end do;
#The results:
seq(res[i],i=1..14);
#Checking results:
seq(q(op(res[i])),i=1..14);

There is no such thing as the particular solution. Any solution is a particular solution. Thus the simplest way to get a particular solution without arbitrary constants is to use dsolve after which one can set the arbitrary constants to whatever one wishes. I shall choose to set both to 0.

restart;
f := diff(u[1](t), t, t)+u[1](t)+(3/4)*a^3*cos(beta+t)+(1/4)*a^3*cos(3*beta+3*t)=0;
dsolve(f);
eval(%,{_C1=0,_C2=0});

This is weird indeed, but has clearly to do with the appearance of a''(t) as a 'known' function.

There is a similar procedure for bvp-problems `dsolve/numeric/bvp/convertsys`.
When you try

`dsolve/numeric/bvp/convertsys`(test, [], x(t), t, y, yp );

you receive the error message

Error, (in dsolve/numeric/bvp/convertsys) specified dependent variables x(t) do not agree with input system {a, x}

A workaround is to spell diff incorrectly as in ddiff(a(t),t,t) or of course just use a notation like app(t).

Once you receive the result from convertsys, you can easily replace those by what they ought to be.
################
The question remains why replacing 'a' by 'p' helps.

Looking at the procedure `DEtools/convertsys` in particular at the loop starting with line 51:
for idiff in diffset do ....
it appears that the order in which the elements of diffset are taken is different in the two cases (a(t) versus p(t)) in Maple 2015 and also in Maple 18, but not in Maple 12 or 15.
I found this confirmed also by
diffset:={diff(a(t), t), diff(x(t), t), diff(a(t), t, t), diff(x(t), t, t)};
convert(diffset,list);
diffset:={diff(p(t), t), diff(x(t), t), diff(p(t), t, t), diff(x(t), t, t)};
convert(diffset,list);
In Maple 2015 (and 18) the a's come first in the first case. In the second case the second derivatives come first.
In Maple 12 and 15 the first derivatives come first and convertsys doesn't work for either use of letter.

You have 'res' in the definition of T0 and T1 inside the procedure Q2. That line should be:

T0,T1:=op(subs(subs(solT),[theta(x),diff(theta(x),x)]));

You may consider simplifying EQT2 by doing:

EQT2:=simplify(EQT2);

That makes EQT2 shorter.

You could do as follows:

Th0:=[seq(cat(Y1,i),i=1..5)](0); #The theta(0) values
PrTh0:=zip(`[]`,L,Th0); #The points in the Pr-theta(0) plane
plot(PrTh0,labels=[Pr,theta(0)]); #Points connected
plot(PrTh0,labels=[Pr,theta(0)],style=point,symbol=solidcircle,symbolsize=20); #Points not connected

#################################
Using an arbitrary number of Pr values:

restart;
  Digits:=15:
  Eq1 := diff(f(eta), eta, eta, eta)+f(eta)*(diff(f(eta), eta, eta))-(diff(f(eta), eta))^2 = 0;
  N := 1;
  blt := 10;
  Eq2 := (diff(theta(eta), eta, eta))/Pr+f(eta)*(diff(theta(eta), eta)) = 0;
  bcs1 := f(0) = 0, (D(f))(0) = 1, (D(f))(blt) = 0;
  bcs2 := (D(theta))(0) = -N*(1+theta(0)), theta(blt) = 0;
# Procedure to compute theta(0) for given value of Pr:

p:=proc(pr) local res;
  res:=dsolve(eval({Eq1, Eq2, bcs1, bcs2}, Pr = pr), numeric, output = listprocedure,maxmesh=1024,initmesh=1024);
  subs(res,theta(eta))(0)
end proc;
 
plot(p,2.5..10);

p:=x^(n-1)+y^g+x^m+5*x^(n-2);
max(op(map2(op,2,indets(p,`^`(identical(x),anything)))))  assuming n>=m+1;

It is sometimes better (in my view) to use linear algebra for linear systems. It provides more insight.

I use the version of your system as corrected by Markiyan.

restart;
sys_ode := diff(ph(t),t) = (1-yc)*pc(t)+yh*prj(t)+urd*prd(t)+ugd*pgd(t)-yc*ph(t), diff(pc(t),t) = yc*ph(t)-(2-yc)*pc(t), diff(pa(t),t) = ya*pc(t)-pf(t), diff(prj(t),t) = yrj*pa(t)+prj(t), diff(prd(t),t) = -urd*prd(t)+yrd*pa(t), diff(pgd(t),t) = -ugd*pgd(t)+ygd*pa(t), diff(pf(t),t) = (1-ygd-yrj-yrd)*pa(t)+(1-yh)*prj(t);

ics := ph(0) = 1, pc(0) = 0, pa(0) = 0, prj(0) = 0, prd(0) = 0, pgd(0) = 0, pf(0) = 0;

VARS:=map2(op,1,lhs~([sys_ode])); #The functions in the order of the derivatives on the left hand sides.
LinearAlgebra:-GenerateMatrix(rhs~([sys_ode]),VARS): A:=%[1];
# Now your system is in vector form:
diff~(Vector(VARS),t)=A.Vector(VARS);
#The eigenvalues of A are the 7 roots of p:
p:=LinearAlgebra:-CharacteristicPolynomial(A,lambda);
# Way too complicated in general.
#So we take a concrete case: All parameters are set to 1.
#But obviously you could just set them to what you want.
params:=indets({sys_ode},name) minus {t} =~1;
eval(p,params);
fsolve(%=0,lambda,complex); #The eigenvalues
#In order to find the solutions we need eigenvalues and eigenvectors:
L,V:=LinearAlgebra:-Eigenvectors(evalf(eval(A,params)));
#Notice that VARS and ics come in the same order:
VARS;
ics;
# The initial value vector:
ICS:=Vector(rhs~([ics]));
# The solution from linear algebra:
sol:=V.LinearAlgebra:-DiagonalMatrix(exp~(t*L)).V^(-1).ICS;
# We know a priori that all the solutions will be real valued functions. We would like to see that they actually are:
evalc(sol);
SOL:=fnormal~(%);
plot(SOL,t=0..5);
# Testing that our solution is actually a solution:
odetest(VARS=~SOL,eval({sys_ode},params));
fnormal(%);
simplify(%);
# If you happen to want the integral of 1-pf(t) from t to infinity you will see that it must be infinite:
int(1-SOL[-1],t=0..infinity);
##############

From your uploaded worksheet it appears that
params:={yc=1/3600,yh=1,yrd=1/180,ygd=1/180,yrj=1/180,urd=0.8,ugd=0.8,ya=1/180};

This won't change the conclusion, however. You may want to set Digits:=15 and plot from t=0..100, but otherwise the code should not need change.

You have list brackets around your expression assigned to f_jp.

It should just be

f_jp := Re(Lf((b+I*li[j]+I*evalf(Pi)*(p/m-1))/Delta));

If you want to keep the big O then you can do:

u1:=taylor(u(x-h,y),h,4);
u2:=taylor(u(x+h,y),h,4);
series(u1+u2,h);

If on the other hand you really don't want the big O anywhere, you can use mtaylor instead:

u1:=mtaylor(u(x-h,y),h,4);
u2:=mtaylor(u(x+h,y),h,4);
u1+u2;

Finally, the shortest possible versions of both:
taylor(u(x-h,y)+u(x+h,y),h,4);
mtaylor(u(x-h,y)+u(x+h,y),h,4);

@Christopher2222 Isn't the sign on U wrong?  (Note: sin instead of cos, se remark below).
If you correct that you could use events to make the ball bounce back at theta = 0:

sol := dsolve({EL, ini},numeric,events=[[theta(t)=0,diff(theta(t),t)=-diff(theta(t),t)]]);

If you add the range argument range=0..10 to dsolve and use refine=1 in odeplot, you get a nicer curve.


There is a similar example in the help page for dsolve/numeric/events.

There is no point in trying to convert the NULL response from your dsolve command to Bessel.

Edited:
Try this:
dsolve({eq1,eq2,eq3}); #General solution: involves BesselJ and BesselI
#We try to take it stepwise:
dsolve({eq1,theta(a)=A,theta(0)=1}); #Doesn't work. Probably because of BesselI(0,k*y) having a singularity at y=0).
#So try using b instead of 0 and taking limit.
res1:=dsolve({eq1,theta(a) = A,theta(b)=1});
TH:=limit(rhs(res1),b=0);
#However, notice that because of the singularity at y=0 of BesselI(0,k*y) merely imposing the condition that theta(0) be finite forces the constant _C2 in the result of
dsolve(eq1);
to be 0.
Thus although TH satisfies theta(a)=A it will not necessarily satisfy theta(0)=1.
In other words you cannot impose more than the one condition at y=a if you ask for a solution that is finite at y=0.

There are several problems:

1. A non-matching parenthesis ) is revealed when copying and pasting your dsys. I removed one of the 5 in ))))).
2. Curly braces {} are used around -1 in what should be z^(-1) in two places.
3. The initial conditions are given at z=0: Not so good because of z^(-1).
4. To solve the system you need to use a DAE method, e.g. method=rkf45_dae.
5. You can have only one condition on u.
6. The DAE-method will differentiate the second ode. The initial conditions must satisfy the second ode (before differentiation).

As an example this works:
restart;
sys := {(1-4*(diff(ln(v(z)), z)))*(diff(u(z), z))+((3/2)*z^(-1)-2*(3* z^(-1) *(diff(ln(v(z)), z))+2*(diff(ln(v(z)),z,z ))))*u(z) = 0, -z*(diff(v(z), z))-v(z)+v(z)^(1/2)*u(z) = 0};
ics:={u(0.1)=1,v(0.1)=1,D(v)(0.1)=0};
res:=dsolve(sys union ics, numeric,method=rkf45_dae);
plots:-odeplot(res,[[z,u(z)],[z,v(z)]],0.1..3);
#################

As for the analytic solution you give I find no requirements needed such as z<0 and z>-4.
sol:= {u(z)=(1+z/4)*exp(z/8), v(z)=exp(z/4)};
odetest(sol,sys);
simplify(%) assuming real;
diff~(sol,z);
convert(%,D);
eval(%,z=0);
eval(sol,z=0);

To find parameter values for which your trial solution is actually a solution, you can use solve/identity as follows.

Begin by defining odesys as you do. Change gamma in your trial solution puta to gamma1 because gamma is Euler's constant in Maple.

Then do:
res := op(odetest(puta, odesys));
SOL:=solve(identity(res=0,r),[a,alpha,beta,gamma1,lambda]);

In Maple 2015 we find 8 different solutions for the parameters.