First a comment:

Your definitions of the functions g1 and g2 only give you the intended result if used with the input 't' (i.e. the letter t).
To see my point just try g1(s); or g1(2);

restart;
x1 := sech(sqrt(a)*t); x2 := -sech(sqrt(a)*t)*tanh(sqrt(a)*t);
Int((Int((eval(x2, t = tau))*sin(Omega*tau), tau = 0 .. u))/(eval(x2, t = u))^2, u = 0 .. t);
g11:=unapply(%,t,Omega,a);
plot3d(rcurry(g11,1), 0 .. 1, 0 .. 3,axes=boxed,labels=[t,Omega,g1],grid=[10,10]);
Int((Int((eval(x2, t = tau))*cos(Omega*tau), tau = 0 .. u))/(eval(x2, t = u))^2, u = 0 .. t);
g22:=unapply(%,t,Omega,a);
plot3d(rcurry(g22,1), 0..1,0..3,axes=boxed,labels=[t,Omega,g2],grid=[10,10]);
#The integrand in your last integral is even:
Ig:=eval(x1*x2^3,a=1)*g11(t,Omega,1)*g22(t,Omega,1);
IntegrationTools:-Change(eval(Ig,t=-t),u=-v);
IntegrationTools:-Change(%,tau=-s);
subs(s=tau,v=u,%);
%-Ig;
simplify(%);
#Thus it suffices to integrate from 0 to infinity and multiply by 2.
# Is the last integral convergent?
The following suggests that it might not be (or that at least it requires a proof of convergence):
Take as an illustration Omega=1.

f:=eval(Ig,Omega=1);
#For large values of t we may do:
eval(f,{tanh=1,sech=2/exp});
simplify(%);
value(%);
g:=combine(expand(%));
map(int,[op(g)],t=0..infinity);

Well the procedure is almost written in the first few lines of your question:

p:= proc (n::posint) option remember;
   if n = 1 then
      1
   else
      p(n-1)+2*n-1   
   end if;
end proc;
p(1);
p(2);
p(6);
#Check (and much easier, but not recursive):
add(2*k-1,k=1..6);
p(45);
p(46);
#You don't need option remember, but in case you use it, here is the content of the remember table so far:
op(4,eval(p));

You didn't square y I see in your picture.
V:=(x,y)->-10/sqrt((x+10/11)^2+y^2)-1/sqrt((x-100/11)^2+y^2);
plot3d(V(x,y),x=-10..15,y=-10..10,view=-10..0,axes=boxed);
plot3d(V(x,y),x=-10..15,y=-10..10,view=-10..0,axes=boxed,style=patchcontour);
plots:-contourplot(V(x,y),x=-10..15,y=-10..10,grid=[50,50],contours=[seq(k,k=-10..-1)]);

Notice first that if you limit arguments to cos to an integer times a name then applyrule could be used:

restart;
rn:=cos(th::name*n::integer)=2*cos(n*th/2)^2-1;
applyrule(rn,cos(alpha));
applyrule(rn,cos(2*alpha));
#However, evalindets might be what you want:
u:=randpoly(x)+t*cos(38+alpha/7);
evalindets(u,specfunc(anything,cos),s->2*cos(op(s)/2)^2-1);
#If you don't want expanded all occurrences of cos, but only those with a specific argument, you could do:
evalindets(u+cos(y),identical(cos(38+alpha/7)),s->2*cos(op(s)/2)^2-1);

Another idea could be to change applyrule so that it only iterates N times. Maybe something like
ap:=proc(_rules,expr,N::posint:=infinity) ..................
......................
......................
for j to N while answer<>i do
..................
end proc:
For the rest see
showstat(applyrule);

Whether that is a good idea or not needs some looking into. In the trivial cases considered it works (also with your original rule r). If a good idea then the case N=0 should perhaps be included as a 'do nothing',i.e.
if N=0 then return expr end if;

Using programmer indexing (= using parentheses instead of brackets) you don't have to declare the size in advance.
restart;
B:=Matrix():
dt := 0.1e-2;
xo := 1; vo := 0; omega2 := 1;
B(1,1):=xo: B(1,2):=vo: B(1,3):=0:
for i from 2 to 1000 do  
   xn := xo+vo*dt;
   vn := vo-omega2*xo*dt;
   B(i,1):=xn;
   B(i,2):=vn;
   B(i,3):=i*dt;
   xo:=xn;
   vo:=vn
end do:

B[1,1..3];
B[1000,1..3];
B[1..10,1..3];
plot(B[1..-1,1..2],labels=[x,v]);
plot(B[1..-1,3],B[1..-1,1],labels=[t,x]);
plot(B[1..-1,[3,1]],labels=[t,x]);

You may take a look at overload:

myproc:=overload(
   [
     proc(a::seq(algebraic),$) option overload;
         `*`(a)
     end proc,
     proc(a::seq({set,list}(algebraic)),$) option overload;
         map(`*`@op,[a])
     end proc
  ]);
myproc(a,b,c,9);
myproc();
myproc({1,3,4},{2,3,4},{5,2,1});

How about this:

S:="xxHhxzahnnnHbbb";
StringTools[SearchAll](["h","H"] , S);

ex:=(-14+9*(2-2*cos(mux))^(1/2))*cos(mu*x)+62+27*(2-2*cos(mux))^(1/2);
(Thanks to Alejandro for typing).
s:=op(indets(ex,`^`(anything,1/2)));
subs(s=freeze(s),ex);
collect(%,freeze(s));
thaw(%);
#Or use frontend:
frontend(collect,[ex,s], [{`+`,`*`},{}]);

You define

G := ZeroMatrix(3*(n-1)+6, 3*(n-1)+6);

thus you cannot change G.

Try
G := Matrix(3*(n-1)+6, 3*(n-1)+6);

instead.

For each given value of b you are looking for a Lambda s.t. there exists a nontrivial solution g to Eq2 with the homogeneous boundary values as given. dsolve easily solves the boundary value problem for Eq1, so that f will be known once b is given.
It is to be expected that if a Lambda is found there will be infinitely many solutions g.
Thus when turning the problem for g into an initial value problem with one parameter taken to be dg2= g''(0) it may be OK to just ask for a solution having e.g. dg2 = 0.1 (as I have done in the worksheet, see below).

For each b it appears that there will be several eigenvalues Lambda (as is demonstrated for b=0.5 in the worksheet).
This complicates making a decent plot of Lambda versus b as there will be several branches (infinitely many?).

But maybe this can be of some help:
(version 4 now)

NonlinearEigenvalu.mw

There is 'eliminate', which you could try. However I found the following slightly neater:

restart;
eq1,eq2:=x = a*(3*cos(t) - cos(3*t)),y = a*(3*sin(t) - sin(3*t));
sys := [eq1,eq2]:
expand(sys);
sys1:=subs(cos(t)=c,sin(t)=s,%);
solve(sys1,{s,c});
eq:=subs(%,s^2+c^2=1);
indets(eq,specfunc(anything,RootOf));
alias(q=op(%));
eq;
normal((lhs-rhs)(%));
EQ:=simplify(numer(%));
EQA:=[allvalues(EQ)];
plots:-implicitplot(eval(EQA,a=1),x=-3..3,y=-4..4,grid=[50,50],gridrefine=1);
plot(eval([rhs(eq1),rhs(eq2),t=-Pi..Pi],a=1));

Well, are you serious? There are over 2000 periods of sin.
You could raise numpoints in effect forcing plot to paint the whole thing red (or whatever color).
restart;
phi:=-223000+2*(118470000-3300*t)^(5/8);
h:=(1/(3590-t))^(1/4)*sin(phi);
evalf((eval(phi,t=0)-eval(phi,t=3600))/2/Pi);
plot(h,t=0..3589,numpoints=5000);

Nothing prevents you from using numpoints=10000 if you are still not satisfied!

You need to write BCS union ICS as they are both sets.
There is a simple typographical error in PDE2. The arguments to theta should be x,t not y,x

It is sometimes easier to solve this by turning it into an initial value problem and by using the parameters option in dsolve. I have often suggested this for systems looking very similar to the one you have.
But here it is:

sol:=dsolve({ode1,ode2,D(f)(0) = 0, f(0) = 0, (D@D)(f)(0) = f2, g(0) = 1, D(g)(0) = g1},
numeric,parameters=[g1,f2],output=listprocedure);
F1,G:=op(subs(sol,[diff(f(x),x),g(x)]));
p1:=proc(g1,f2) sol(parameters=[g1,f2]); F1(6) end proc;
p2:=proc(g1,f2) sol(parameters=[g1,f2]); G(6) end proc;
p1(0.4,1);
p2(0.4,1);
fsolve([p1,p2],[0.4,1]);
sol(parameters=%);
plots:-odeplot(sol, [[x,diff(f(x),x)],[x, g(x)]], 0 .. 6, color = [blue,black]);
#Check:
F1(6);
G(6);
#And now that we have a solution we could use that as an 'approximate sloution':
dsn := dsolve(sys, numeric,approxsoln=sol);
plots:-odeplot(dsn, [[x,diff(f(x),x)],[x, g(x)]], 0 .. 6, color = [blue,black]);

If I understand you correctly then the following may be useful, and you don't have a pde, but an ode.
I'm nowhere near solving the boundary value problem (one of the conditions isn't clear to me as presented and the other is given at x = -infinity).
There are enough difficulties with the initial value problem.

restart;  
ode:=U*diff(F(s(x)),x)= -k*sigma*cos(theta)*diff(F(s(x))*kro(s(x))*diff(F(s(x)),x), x)/muo/sqrt(k/phi);
krw:=s->A*(s-swi)^(aw);
#Defining the known functions:
kro:=s->C*(1-sor-s)^(ao);
F:=s->krw(s)/muw/(krw(s)/muw+kro(s)/muo);
J:=s->alpha*(((s-swi)/(1-sor-swi))^(-1/m)-1)^(1-m);
#The next step is not really necessary, but is done to gain understanding of the ode.
ode1:=factor(isolate(ode,diff(s(x),x,x)));
#The numerical constants:
krowi := 1; krwor := .1; aw := 1.16; ao := 3.69; alpha := 0.6e-1; U := 2.9374109836-10^(-7); L := 0.552566e-1; m := .5; sor := .3; swi := .25; k := 10^(-14); muw := 0.1e-2; muo := 0.5e-2; sigma := 0.72e-1; phi := .3; theta := 0;
C := krowi*(1-sor-swi)^(-ao);
A := krwor*(1-sor-swi)^(-aw);
#Taking a look at ode1:
ode1;
ode2:=factor(ode1);
#In order that all the fractional powers involved return real values you need s to satisfy
# swi < s < 1- sor
#Now an attempt:
sol:=dsolve({ode2,s(0)=.6,D(s)(0)=1e6},numeric,stiff=true);
sol(1);
plots:-odeplot(sol,[x,s(x)],0..2,thickness=2);

#This is as far as I got.