After correcting the syntax in your initial conditions you still have a problem.

You have two equations determining diff(psi(t),t,t), eq1 and eq3. Together with eq2 they seem to be inconsistent.

restart;
eq1:=(l^2+r^2/2)*(sin(theta(t))^2*diff(diff(psi(t),t),t)+2*diff(psi(t),t)*diff(theta(t),t)*cos(theta(t)))-r^2*(cos(theta(t))*diff(psi(t),t)+omega(t))*diff(theta(t),t)=0;
eq2:=(l^2+r^2/2)*(diff(diff(theta(t),t),t)-(diff(psi(t),t))^2*sin(theta(t))*cos(theta(t)))+r^2*(diff(psi(t),t)*cos(theta(t))+omega(t))*diff(psi(t),t)*sin(theta(t))=g*l*sin(theta(t));
eq3:=diff(diff(psi(t),t)*cos(theta(t))+omega,t)=0;
l:=0.05;r:=0.063/2;g:=0.81;omega:=evalf(2*Pi*10000/3600);
ic1:=theta(0)=evalf(9*Pi/180),psi(0)=0,D(theta)(0)=0,D(psi)(0)=0;

eq1;eq2;eq3;

#This doesn't succeed:
Sol:=dsolve({eq1,eq2,eq3} union {ic1},{theta(t),psi(t)},numeric);
Maybe there is a good reason. From eq1,eq2, and eq3 you can find two odes, one of which is second order in theta (L2) and the other (L13) is first order in psi. These two can be solved: 
L1:=op(expand(solve(eq1,{diff(psi(t), t, t)})));
L3:=op(solve(eq3,{diff(psi(t), t, t)}));
L2:=op(solve(eq2,{diff(theta(t), t, t)}));
eval(L1,L3);
(rhs-lhs)(%)=0;
L13:=op(solve(%,{diff(psi(t), t)}));
ic1;
ic2:=ic1[1..3];
Sol:=dsolve({L13,L2} union {ic2},numeric,output=listprocedure);
#However, it appears that the second derivative of psi found from Sol by fdiff and the same found from L3 using Sol don't agree:
TH,dTH,PSI:=op(subs(Sol,[theta(t),diff(theta(t),t),psi(t)]));
fdiff(PSI,[1,1],[1.2345]);
psi2:=eval(rhs(L3),L13);
eval(psi2,{theta(t)=TH(t),diff(theta(t),t)=dTH(t)});
eval(%,t=1.2345);

Your f is not defined as a function, so should not be used as such.

Replace f(x) by f:

s := solve(y = f, x);

#Maple reports 9 solutions:
nops([s]);

#one of which contains no obviously imaginary parts:

remove(has,[s],I);

Since DEtools is not a module you cannot use 'uses' and you should not use 'with' within procedures.

type(DEtools,`module`);
type(DEtools,table);

DEtools is a table (used for packages before modules were introduced in Maple).
Use the long form DEtools[DEplot] inside procedures.

Example.

eq:=diff(y(t),t)=y(t);

p:=proc(t0::realcons,y0::realcons,T::realcons) global eq,y,t;    
    DEtools[DEplot](eq,y(t),t=t0..T,[y(t0)=y0],_rest)
end proc;

p(0,1,2,linecolor=blue);

By experience I know that one has to be extremely careful when using an indexed name like 'a[12]' and also its unindexed form 'a'. However, when they are not assigned values, it should not be as dangerous.

Nevertheless, to avoid speculating about that I replaced the bare form 'a' with 'b'.

I ran the following, and found solutions in terms of RootOf's, which is not surprising. But I only found valid solutions when first converting tanh to exp.

restart;
sys:=[-xt[1,1]+b*tanh(xt[1,1])+a[12]*tanh(xt[2,1])=0,-xt[2,1]+a[21]*tanh(xt[1,1])+b*tanh(xt[2,1])=0];
solve(sys,[xt[1,1],xt[2,1]]);
allvalues(%);
#Solutions false:
eval(sys,{xt[1,1]=ln(I),xt[2,1]=ln(I)});
convert(sys,exp);
res:=solve(%,[xt[1,1],xt[2,1]]);
indets(res,name);
#Testing the solution for concrete randomly chosen parameters in the interval (0,10)
#and looking for nontrivial solutions only.
rnd:=10*RandomTools:-MersenneTwister:-GenerateFloat:
RES0:=[[xt[1,1]=0.,xt[2,1]=0.]]: RES:=RES0:
for k from 1 to 100 while RES=RES0 do
  val:=[b=rnd(),a[12]=rnd(),a[21]=rnd()];
  RES:=simplify(fnormal(evalf(eval(res,{val[1],val[2],val[3]}))))
end do:
RES; val;
k;
eval(sys,{op([1,1..2],RES),op(val)});

To write g_{i,j} is certainly not allowed in 1d-input, which I always use.

But if you do

p:=(k,j)->if k=0 then max(k-j*S,0) else g(k,j) end if;

after having defined the function g in the proper manner (*), then p(1,4) should execute as expected.

(I changed K to k assuming that was what you meant).

(*) Example. g:= (i,j) -> i^2-j*i;

The derivative exists and is 1/2.

A:=Int(cos(1/sqrt(u))^2,u=0..x);

IntegrationTools:-Change(A,u=t^2);

B:=combine(convert(%,radical));

#We need only consider

C:=Int(t*cos(2/t), t = 0 .. sqrt(x));

IntegrationTools:-Change(C,t=sqrt(x)/u,[u]) assuming x>0;

#Now we shall show that
E:=Int(cos(n*u)/u^3, u = 1 .. infinity);

#tends to zero as n tends to infinity.

#We split the integral in two, where M is to be chosen later
IntegrationTools:-Split(E,M);

# and use an integration by parts on the first term:
applyop(IntegrationTools:-Parts,1,%,u^(-3));

The remaining details follow by a standard argument.

There are a few errors in your statement of the problem. sin*t should be sin(t), 'pi' should be 'Pi', u(0,t) should not be assigned to. 

These errors are corrected below.

PDE:=diff(u(x,t),t)+2*diff(u(x,t),x)=0;
bc:=u(0,t)=piecewise(0<t and t<=2*Pi, sin(t),0); 
pdsolve({PDE,bc});
U:=rhs(%);
plot(eval(U,t=0),x=-5*Pi..Pi,thickness=2,caption="The solution for t = 0");
plots:-animate(plot,[U,x=-5*Pi..7*Pi,thickness=2],t=0..3*Pi);

Using the polar form for t, t = r*exp(I*v), we immediately see that v = Pi/5 + p*2*Pi/5, where p = 0, 1, 2, 3, 4.

Then

restart;
eq:=(-t)^5=abs((t^2+t+1)/(-t));
evalc(eval(eq,t=r*exp(I*(Pi/5+p*2*Pi/5)))):
eq3:=combine(%) assuming p::integer,r>0;
#Trying p=0:
solve({eval(eq3,p=0),r>0},r);
evalf(%);
#Now trying all p's:
[seq(solve({eq3,r>0},r),p=0..4)]:
L:=evalf(%);
nops(L);
L2:=subs~(L,r);
#Maple found nothing for p = 2, which corresponds to t = -1.
#We consider p = 2 separately:
solve({eval(eq3,p=2),r>0},r);
factor(eval(eq3,p=2));
simplify(%) assuming r>0;
solve({%,r>0},r);
res:=zip(`*`,L2,[seq(exp(I*(Pi/5+p*2*Pi/5)),p=[0,1,3,4])]);
plots:-complexplot([-1,op(res)],style=point,symbolsize=20);

It seems from Maple that the only solution is -1.

I did the following:

eq:=(-t)^5=abs((t^2+t+1)/(-t));
res:=solve(eq,t);
op(0,res);
evalf(res);
#None of the following are satisfied for obvious reasons!
map(evalf,op(1,res));

#The seven values:
evalf(op(2,res))

The answer from solve is certainly not perfect. It comes in the form of a piecewise expression.

Maple finds 7 t-values (6 RootOf's and -1). Maple stalls on the condition that the fifth powers of the RootOf's be negative, which they aren't: They are not even real.

Apart from the stalling, if Maple is right, then the only solution is -1.

The problem is that y is used in two different ways. y is first implicitly defined as a table by the assignment y[1]:=0.5:

Later it is used as just a name in

s:=solve(y=g[i](x),x):
h[i]:=unapply(s[1],y);

If you change 'y' to e.g. 'yy' the loop runs.

The use of a table name as a mere name is likely to cause strange problems. The resulting error messages are usually very opaque.

If you could provide us with the actual system and initial conditions then we wouldn't have to make guesses.

Using that the two mixed partials w.r.t. x and y are equal you can write down a condition for solving the system:

restart;
pde1 := diff(f(x,y),x) + A(x) * ( p11* f(x,y) + p12 + p13*x )
= p14 * B(y)/B(x) * diff(f(x,y),x);

pde2 := diff(f(x,y),y) + A(y) * ( p21* f(x,y) + p22 + p23*y )
= p24 * B(x)/B(y) * diff(f(x,y),y);

diff(isolate(pde1,diff(f(x,y),x)),y);
diff(isolate(pde2,diff(f(x,y),y)),x);
factor(rhs(%%-%));
Nm:=numer(%):
#Your first special case:
collect(eval(Nm,{p11=0,p21=0}),diff,factor); #This has to be zero
#Your second special case:
factor(eval(Nm,{p11=0,p21=0,p24=0}));#This has to be zero

In both cases you get requirements on the coefficients for solutions to exist.

In the latter special case the requirement is that

A(x)*B(y)^2*B(x)*(diff(B(y), y))*p14*(p12+p13*x)=0 for all (x,y) in some open non-empty set.

Unless p14 = 0 or p12 = p13 = 0 then the only interesting case would be B = constant.

All parameters must have values, and you need one more initial condition.

I have added D(r)(0)=0.

restart;
G := 1;
M := 1;
pT := 2;

OrbitDEs := diff(r(tau), [tau$2]) = -G*M/r(tau)^2+L^2/r(tau)^3-3*G*M*L^2/r(tau)^4, diff(phi(tau), tau) = L/r(tau)^2, diff(t(tau), tau) = sqrt(r(tau)/(r(tau)-2*G*M)+(diff(r(tau), tau))^2*r(tau)^2/(r(tau)-2*G*M)^2+r(tau)^3*(diff(phi(tau), tau))^2/(r(tau)-2*G*M));

OrbitICs := r(0) = r0, phi(0) = 0, t(0) = 0, D(r)(0)=0;
#Arbitrarily using L = 1, r0=3, and pT= 7.
#You don't have to provide the range argument range = 0..pT.


eqs:=eval({OrbitDEs, OrbitICs},{L=1,r0=3});
OrbitSolution := dsolve(eqs, numeric, {r(tau), phi(tau), t(tau)}, range = 0 .. 7);

plots:-odeplot(OrbitSolution,[tau,phi(tau)],0..3.9);

The error message you reported said to take a look at dsolve,numeric,parameters. That is a good idea, if you intend to change the values of the parameters.

There are two complex solutions to z^2 = a, where a is a given complex number, i.e. two square roots if you will.

Maple's sqrt will find the square root having the least principal argument in absolute value, i.e. the principal square root. (The same is true for other roots).

One square root is sqrt(a) the other is -sqrt(a).

The equation y(lambda,beta) = 0 can be turned into an ode problem

For simplicity:
y:=subs(`&beta;b`=b,y);

eq:=diff(b(lambda),lambda)=-eval(diff(y,lambda)/diff(y,b),b=b(lambda)):

It appears that y is purely imaginary for b between 1 and sqrt(2.1025). For b<1 or b > sqrt(2.1025) and b<sqrt(12.25) y is real.

In the latter two regions the ode can be solved numerically. The initial condition can be found using fsolve with the proper initial point.