You could try using the result from dsolve/numeric/bvp to find the initial values at one of the two endpoints, in your case they are -0.5 and 0. Indeed you could use any point inside the interval instead. Then use these as initial conditions in dsolve/numeric.
I have tried this in an extremely simple case:

ode:=diff(y(x),x,x)+y(x)=0;
res:=dsolve({ode,y(0)=0,y(Pi/2)=1},numeric); #This is the bvp solution
plots:-odeplot(res,[x,y(x)-sin(x)],0..Pi/2); #Just checking how good it is.
ics:=eval(convert(res(0),D),x=0)[2..-1]; #Will be taken as initial conditions
res_ics:=dsolve({ode,op(ics)},numeric); #NOT using dsolve/numeric/bvp
plots:-odeplot(res_ics,[x,y(x)-sin(x)],-2*Pi..3*Pi); #Not too bad

I find it a bad idea to use lower case l, since it looks quite like 1.

You have

(D(f))(0) = l+b*(D@@2)(f)(0);

where l := [1, 2, 3], which clearly is a problem.

I don't know why, but dsolve/numeric/bvp seems to run into a problem somewhere below 0.40.
I don't know if it is a genuine problem or a technicality, but you can try this:

res:=dsolve(eval([sys_ode, ics], n = .3787*c),numeric, output=operator,
       continuation = c,maxmesh=2048);
This continuation process starts with c=0 and ends with c=1, i.e. starts with n=0 and ends with n=0.3787 as I have it above. Surely it works by using the previous result as an approximate solution for the present.
Incidentally, you could make that process yourself in a loop by using the optional argument approxsoln=res, where res is the solution found previously. Warning: Don't use output=operator for that process as that currently doesn't work. Just change to output=listprocedure (or the default, procedurelist).

Here is a solution using dsolve/numeric/bvp on the two intervals -0.5..0 and 0..5.
Two parameters are introduced. The parameters are determined by fsolve by requiring smoothness.
You can try the other values in your list a yourself. I show the result for a[1], which is referred to here as A.
The procedure p is primarily for use in finding the parameters, but will also serve to present the result, when output is set at any other name than 'number' (the default).

restart;
sys := diff(f(eta), eta, eta, eta)+f(eta)*(diff(f(eta), eta, eta))+1-(diff(f(eta), eta))^2 = 0, diff(theta(eta), eta, eta)+.72*(diff(theta(eta), eta))+A*(diff(f(eta), eta))+theta(eta) = 0;
bcs := f(-.5) = 0, (D(f))(0) = 1+((D@@2)(f))(0), (D(f))(5) = 0, theta(-.5) = 1+(D(theta))(0), theta(5) = 0;
A:=-.5:
p:=proc(f1,th1,{output::name:='number'}) option remember; local res1,fvals,thvals,res2;
   res1:=dsolve({sys,f(-.5) = 0,D(f)(0) = f1, (D@@2)(f)(0)=f1-1,theta(-.5) = 1+th1,D(theta)(0)=th1},numeric,:-output=listprocedure);
   fvals:=subs(res1,[seq(diff(f(eta),[eta$i]),i=0..2)])(0);
   thvals:=subs(res1,[seq(diff(theta(eta),[eta$i]),i=0..1)])(0);
   ##
   res2:=dsolve({sys,f(0)=fvals[1],D(f)(0) = fvals[2],theta(0)=thvals[1],D(f)(5) = 0,theta(5) = 0},numeric,:-output=listprocedure);
   if output='number' then
      [fvals[3]-subs(res2,diff(f(eta),eta$2))(0),thvals[2]-subs(res2,diff(theta(eta),eta))(0)]
   else
       res1,res2
   end if
end proc;
p1:=proc(f1,th1) p(_passed)[1] end proc;
p2:=proc(f1,th1) p(_passed)[2] end proc;

p(.3,-.2); #Just testing the procedure p
par:=fsolve([p1,p2],[.3,-.2]);
res1,res2:=p(op(par),output=xxx);
plots:-display(plots:-odeplot(res1,[[eta,f(eta)],[eta,theta(eta)]]),plots:-odeplot(res2,[[eta,f(eta)],[eta,theta(eta)]]));
plots:-display(plots:-odeplot(res1,[[eta,diff(f(eta),eta)],[eta,diff(theta(eta),eta)]]),plots:-odeplot(res2,[[eta,diff(f(eta),eta)],[eta,diff(theta(eta),eta)]]));
plots:-display(plots:-odeplot(res1,[[eta,diff(f(eta),eta,eta)]]),plots:-odeplot(res2,[[eta,diff(f(eta),eta,eta)]]));

I assume that the norm you will use is the 2-norm over -3..3.
If so do:

soln1 := [seq(dsolve(schro1[i] union Ic[i], {psi(x)}, numeric,output=listprocedure), i = 1 .. nops(E))]; #output=istprocedure

for i from 1 to nops(E) do N[i]:=evalf(Int(subs(soln1[i],psi(x))(x)^2,x=-3..3,epsilon=1e-8))^(1/2) end do;

Then plot

display(seq(odeplot(soln1[i], [x, psi(x)/N[i]], -3 .. 3, color = [red, blue, green][i]), i = 1 .. nops(E)));

Use add instead of sum. In Maple 2015 and 2016 you can tell add which increment you want, here 2:

add((2*q*cos(2* i*x)*(-1)^(i)*(-1)^((2*i-1)))/(i*Pi),i=1...35,2);

If your Maple version is 18 or earlier you can use seq combined with convert/+:

seq((2*q*cos(2* i*x)*(-1)^(i)*(-1)^((2*i-1)))/(i*Pi),i=1...35,2);
convert([%],`+`);

Allowing ourselves to multiply by u' we get:

ode:=diff(u(x),x,x)=a0+a1*u(x)+a2*u(x)^2;
diff(u(x),x)*ode;
map(int,%,x)+(0=C);
solve(%,{diff(u(x),x)});

Using
Digits:=15;

and your system SYS together with inspiration gathered as described in my last comment about the problem, I came up with an approximate solution resembling the graphs in the comment.
Also I used abserr=1e-3.

The following code produced the graphs below. A warning about raising Digits to approximately 16 came up. Then I tried raising Digits to 20 and got a warning about raising Digits to 21, etc. ; a wild goose chase that I have encountered (too) often with dsolve/numeric/bvp. I gave up chasing the goose at Digits:=100.

solA:=dsolve(SYS,numeric,method=bvp[middefer],
   approxsoln=[F(eta)=10^(-3)*eta^2,K(eta)=10^(-6)*tanh(eta),Omega(eta)=2*10^(-9)*tanh(eta),theta(eta)=1-eta], abserr=1e-3);

#########################################################
## Trying a revised set of boundary conditions: Replacing (D@@2)(F)(1)=0 with F(1)=1.
## That I set F(1) = 1 is just taken as an example of setting F(1) at some positive value at 1.
##
sys,bcs:=selectremove(has,SYS,diff);
bcsF1:=subs(( (D@@2)(F)(1)=0 ) = (F(1)=1), bcs); #New boundary conditions!!!
solF1a:=dsolve(sys union bcsF1,numeric,method=bvp[middefer],approxsoln=[F(eta)=eta^2,K(eta)=2*tanh(eta),Omega(eta)=10^(-16)*tanh(eta),theta(eta)=1-eta]);
##
plots:-odeplot(solF1a,[seq([eta,fu(eta)],fu=[F,K,10^16*Omega,theta])],0..1,legend=[F,K,10^16*Omega,theta],color=[red,blue,green,black]);

Using dsolve you can do:

res:=dsolve(omega*x+sigma*y(x)+(omega*x-sigma*y(x))*diff(y(x),x)=0,implicit);
f:=eval(lhs(res),{_C1=0,y(x)=y});
factor(diff(f,x));
factor(diff(f,y));
 from which an integrating factor is clearly seen to be -1/(omega*x^2+omega*x*y+sigma*x*y-sigma*y^2).

If you make the modifications appearing here it works:

interim:=t->evalf(Int(x->r1_lab_deriv(x)[1],0..t,epsilon=1e-6));
interim(5);

Notice:
1. Int instead of int to use numerical integration.
2. You were not integrating a procedure as I would understand it:
  This would be an example of integrating a procedure numerically:
  evalf(Int(sin, 0..Pi));
  thus I have as first input to Int the procedure x->r1_lab_deriv(x)[1], and the second is 0..t (no x).
3. I have added the optional argument concerning accuracy epsilon=1e-6 to get an answer fast.
    

Here is the solution briefly described in my comment:

restart;
eq:=f(s) = 72*s*(-s^2+1)/(Pi*(s^2+1)^4)+Int(f(t)*ln(abs((s+t)/(s-t))), t = 0 .. 1)/(I*Pi);
n:=10:
basis:=[seq(s^i,i=0..n)]; #Choice of basis functions
Y:=unapply(add(c[i]*basis[i],i=1..n+1),s); #To be used for f(s)
eqs:=eval(eq,f=~Y);
## We shall evaluate eqs at n+1 values of s, since we have n+1 unknowns c[i].
#The choice made here is to use the zeros of ChebyshevT(n+1,x) on -1..1, but mapped onto 0..1 as seen in this #simple procedure:
####
ChebAbscissae:=proc(r::range(realcons),n::posint) local a,b,T;
  a,b:=op(evalf(r));
  T:=[seq(cos((2*i+1)/(n+1)*Pi/2),i=n..0,-1)];
  evalf(T*(b-a)/2+~(a+b)/2)
end proc;
####
T:=ChebAbscissae(0..1,10); #We need these collocation points
nops(%);
#### The following 3 lines just serves as a comment on the zeros on -1..1  #################
fsolve(simplify(ChebyshevT(11,x))=0,x);
nops([%]);
evalf([seq(cos((2*i+1)/(n+1)*Pi/2),i=n..0,-1)]);
################################################################
eqs2:={seq(eval(eqs,s=s0),s0=T)}:
eqs3:=subs(Ln=ln,IntegrationTools:-Expand(subs(ln=Ln,eqs2))): #Prevent unwanted expansion of ln
eqs4:=evalf(eqs3); #Numerical integration of the integrals
res1:=fsolve(eqs4); #Solving for the coefficients c[i]
res:=f(s)=eval(Y(s),res1);
plots:-complexplot(rhs(res),s=0..1);

#You may then continue as in my comment with
test:=eval((lhs-rhs)(eq),f=unapply(rhs(res),s));
numapprox:-infnorm(abs(test),s=0..1);
plot(abs(test),s=0..1);

Have you tried
simplify(u,size);

on your expression?

Using the code in Solve1 in the link (and as given in your worksheet):

A:=a^%T;
B:=a;
C:=-q;
X:=Matrix(rtable_dims(A), symbol= x);
eqs := convert(A . X+X . B=~ C, set);
## Output:
eqs := {0 = -500, 0 = -350, 0 = 0, x[1, 1] = 0, x[1, 2] = 0, x[1, 3] = 0, x[2, 1] = 0, x[2, 2] = 0, x[2, 3] = 0, x[3, 1] = 0, x[3, 2] = 0, x[3, 3] = 0, x[1, 4]+x[4, 1] = 0, x[1, 5]+x[4, 2] = 0, x[1, 6]+x[4, 3] = 0, x[2, 4]+x[5, 1] = 0, x[2, 5]+x[5, 2] = -20, x[2, 6]+x[5, 3] = 0, x[3, 4]+x[6, 1] = 0, x[3, 5]+x[6, 2] = 0, x[3, 6]+x[6, 3] = -20}

i.e. there are two equations 0=-500 and 0=-350 that are never satisfied no matter what X is.

To get another view, just look at the output from:
A.X+X.B=C;

Your problem is that you are using a1(I1,E2,I2) and a3(E1,E2,I2), i.e. you have different dependencies in the same system. I suppose you could just expand both to e.g. a1(I1,I2,E1,E2) or som fixed order.
## In fact that works fine:
sys4:=subs(a1(I1,E2,I2)=a1(I1,I2,E1,E2),a3(E1,E2,I2)=a3(I1,I2,E1,E2),sys2);
pdsolve(sys4);
  

I first tried a different approach.

I tried solving your system sys2 for the derivatives:

res:=solve(sys2,{diff(a1(I1, E2, I2), E2), diff(a1(I1, E2, I2), I1), diff(a1(I1, E2, I2), I2), diff(a3(E1, E2, I2), E1), diff(a3(E1, E2, I2), E2), diff(a3(E1, E2, I2), I2)});
nops([res]);

and found 3 sets of results. Then I tried pdsolve on each res[i] for i=1..3.
The inconsistency warning came in all 3 cases for the reason mentioned in the beginning.
For simplicity right now I shall just consider res[1], which I shall call sys:

restart;
sys:={diff(a1(I1, E2, I2), E2) = -sqrt(-E1*I1-E2*I2)*I2/I1, diff(a1(I1, E2, I2), I1) = -(1/3)*sqrt(-E1*I1-E2*I2)*(E1*I1-2*E2*I2)/I1^2, diff(a1(I1, E2, I2), I2) = E2*(E1*I1+E2*I2)/(I1*sqrt(-E1*I1-E2*I2)), diff(a3(E1, E2, I2), E1) = -(1/3)*sqrt(-E1*I1-E2*I2)*(E1*I1-2*E2*I2)/E1^2, diff(a3(E1, E2, I2), E2) = I2*(E1*I1+E2*I2)/(E1*sqrt(-E1*I1-E2*I2)), diff(a3(E1, E2, I2), I2) = -sqrt(-E1*I1-E2*I2)*E2/E1};
pdsolve(sys); #The warning
Now separate the equations in a1 and a3:
K3:=select(has,sys,a3);
K1:=select(has,sys,a1);
pdsolve(K3);
pdsolve(K1);
## Both are successful.

It is not upper case L, but lower. So it is
local aaa,ii;