1. You forgot to assign to a.
2. Use unapply in defining psi_5 .. psi_8:
psi_5:=unapply(psi_1(x,y)*chi(x,y),x,y);
psi_6:=unapply(psi_2(x,y)*chi(x,y),x,y);
psi_7:=unapply(psi_3(x,y)*chi(x,y),x,y);
psi_8:=unapply(psi_4(x,y)*chi(x,y),x,y);

3. Wrap the definitions of the relevant matrices in evalf so that the integrals get computed.
To see if E is a numerical matrix you can do:
interface(rtablesize=infinity);
E;

4. The computation should be fast now, but notice that the datatype of E is anything:
op(E);
To change that you could do E:=Matrix(E,datatype=float);

The results of doing what I described in my last comment above made me think that an approximate solution for the system in the scaled variables y and psi would have
psi[1](t) = N-t
psi[3](t) = 0
y[1](t) = 1
y[2](t) = 1

with psi[2] and y[3] to be determined.
Inserting the 4 "known" functions into the system one finds very simple differential equations for psi[2] and y[3].
Using the conditions psi[2](N) = 0 and y[3](0)=1 dsolve (without 'numeric') finds easily that
psi[3](t) = 0
and
y[3](t) = 2.1891-1.1891*exp(-.67*t).

OK, so how good is this approximate solution?
As far as satisfying the system of odes measured by the difference between rhs and lhs of the equations the answer is that the largest difference when N=100 is 0.0018.

MaplePrimes13-12-07b.mw

z is a global variable in the first and second session. The assignment b:=2*z is saved.
In the third session z is declared local and assigned to k. The global version of b is read. w is assigned to b, which was assigned to the global version of z. Thus w will not be affected by assignments to the local z.

Try changing the declarations to
local w; global z;

The difference is exhibited here:

restart;
ilog[3](3^22-1); #output 22
restart;
Digits:=50:
ilog[3](3^22-1); #output 21

Try

showstat(ilog);

You will see that in the indexed case, using ilog[b], where the index b is 10 or 2, then the builtin functions ilog10 or ilog2 are used, respectively.
In all other indexed cases we have b1 := evalf(b); immediately introducing results possibly depending on the setting of Digits.

V:=1-2*M*exp(-2*l*r)+q^2*exp(-4*l*r);
plot(eval(V,{M=.1,l=.4,q=.7}),r=0..3,filled,color=blue);

Change
part:= int(list[i]*f*(x+1),x=-1..1)*list[i]:
into
part:= int(list[i](x)*f*(x+1),x=-1..1)*list[i](x):

since 'list' is a list of procedures (functions).

Then

F:=functions(3);
returns 3 functions and you can do e.g. F[3](5);

Alternatively (and better, I think) wait with unapply until the very end:

functions:= proc(n)
local L, list, p, f, sum, i, part, g, normg, x:
L:=1/sqrt(2):
list:=[L]:

for p from 2 to n do
  f := x**(p-1):
  sum := 0:

  for i from 1 to (p-1) do
    part:= int(list[i]*f*(x+1),x=-1..1)*list[i]:
    sum := sum + part:
  end do:

  g:= f-sum:
  normg:= sqrt(int(g^2*(x+1),x = -1..1)):
  L:=g/normg:
  list := [op(list), L]:
end do:
unapply~(list,x)
end proc:

You can now do as before.

Let u be your expression. Then try:

evalindets(u,radical,factor);
simplify(%) assuming positive;
q:=denom(%)/lambda;
simplify(algsubs(q =w,%%));
expand(%);
eval(%,w=q);
map(simplify,%) ;
combine(%) assuming positive;

I'm not following you. I cannot see that f1 and f2 are solutions:

odetest(y(t)=f(t),ode2); #Yes, y(t)=f(t) is a solution
odetest(y(t)=f(t-t1),ode2);
eval(%,{t=0,t1=1,s=1});
evalf(%);
#Not zero meaning that y(t)=f(t-t1) is not a solution of ode2.

However, reduction of order gives you a solution linearly independent of f(t) easily:

simplify(eval(ode2,y(t)=f(t)*v(t)));
eval(%,diff(v(t),t)=w(t));
dsolve(%);
combine(int(eval(rhs(%),_C1=1),t)); #Finding v(t)
f(t)*%;
convert(%,trigh);
f2:=simplify(%); #The second solution
odetest(y(t)=f2,ode2); #Indeed
#There is really no reason to check the wronskian since v(t) is not constant in t.
#But if we do:
VectorCalculus:-Wronskian([f(t),f2],t,determinant)[2];
convert(%,exp);
#we find (1/64)*s^2.
# s= 0 is a special case and not very interesting:
eval(ode2,s=0);

You cannot have infinity as the right "endpoint". Choose a large number in this case I chose 10.
I assume that your last boundary condition is meant to mean f''(infinity) = 0. With infinity = N that would be
(D@@2)(f)(N) =0 or (D@D)(f)(N)=0.
Try this then.
N := 10;
bcs := (D(f))(0) = 1, f(0) = 0, (D(f))(N) = 0, ((D@@2)(f))(N) = 0;
sol := dsolve({bcs, eq1}, numeric, method = bvp[midrich], maxmesh = 1024);
odeplot(sol,[eta,f(eta)],0..N);

u:=BesselJ(1,rho*exp(I*Pi/4));
diff(Re(u),rho); # Notice abs(1, ...)
eval(%,rho=0.5);
evalf(%);
# abs(1,w) is the derivative of abs(w)
#try
diff(abs(w(rho)),rho);
# What you are after may be what you get by postponing using Re:
u1:=diff(u,rho);
eval(u1,rho=0.5);
evalf(%);
Re(%);

There is probably no closed form available in the general case.
In two very special cases, yes:

dsolve(eval(ODE[4],a=0));
dsolve(eval(ODE[4],epsilon[2]=0));

Otherwise, solve numerically by maybe using the parameters option in dsolve/numeric as illustrated here.
You need initial conditions. I have used y = 0, but included the values as parameters Y0 and Y1:

res:=dsolve({ODE[4],Y(0)=Y0,D(Y)(0)=Y1},numeric,parameters=[Y0,Y1,a, beta, omega, epsilon[1], epsilon[2], mu[0]]);
#Now set values for the parameters. With this many parameters it is difficult to remember the order, but instead of giving numbers in the correct order you can give equations (in any order) as I have done here:

res(parameters=[Y0=1,Y1=0,a=.1, beta=.2, omega=.3, epsilon[1]=.4, epsilon[2]=.5, mu[0]=-6]);
#This would have been just as good:
res(parameters=[a=.1, beta=.2, omega=.3, epsilon[1]=.4, epsilon[2]=.5, mu[0]=-6,Y0=1,Y1=0]);
#Now plotting
plots:-odeplot(res,[y,Y(y)],0..25);

Q:=Matrix([[a+b,c+d],[e+f,g+h]]);
M:=Matrix([[A,B],[C,F]]);
assign~(M=~Q);
M;
A,B,C,F;

restart;
ode0:=diff(f(x),x,x)/((1+diff(f(x),x)^2)^(3/2)) = e*(x*sin(alpha)+f(x)*cos(alpha))+P;
ics:=D(f)(0)=0, f(0)=1;
params:={e=1/10,alpha=Pi/3};
ode:=eval(ode0,params);
res:=dsolve({ode,ics},numeric,output=listprocedure,parameters=[P]);
F,F1:=op(subs(res,[f(x),diff(f(x),x)]));
K:=1/sqrt(3); #So f'(a) = -1/sqrt(3)
#The procedure p surves dual purposes.
#The default is the primary: Finding f(a) and f'(a)+K.
#The secondary purpose is for use in exploring (plotting).
p:=proc(a,P,output::{procedure,list(procedure)}:='primary') local fa,f1a,i; global p1,p2;
   res(parameters=[P]);
   if output<>'primary' then return output end if;
   p1(_passed):=F(a);
   p2(_passed):=F1(a)+evalf(K);
   if type(procname,indexed) then cat(op(0,procname),op(procname))(_passed) else [p1,p2](_passed) end if
end proc;
p1:=proc(a,P) p[1](_passed) end proc;
p2:=proc(a,P) p[2](_passed) end proc;
#Hinting at exploration
plot(p(.9,-.7,[F,F1]),0..1.5);
aa:=1.9:
plots:-animate(plot,['p(aa,P,[F,F1])',0..aa,-K..1],P=-.7..-.06);
#Solving
sol:=fsolve([p1,p2],[1.9,-.6]);
plots:-odeplot(res,[[x,f(x)],[x,diff(f(x),x)]],0..sol[1]);
F(sol[1]);
F1(sol[1])+evalf(K);

I suggest that you give us the actual ode. In the meantime here is a very simple example, where it is easy to find solutions by hand:
diff(f(t),t,t)+omega^2*f(t)=0.
However, I solve numerically.

restart;
ode:=diff(f(t),t,t)+omega^2*f(t)=0;
res:=dsolve({ode,f(0)=f0,D(f)(0)=0},numeric,output=listprocedure,parameters=[f0,omega]);
F,F1:=op(subs(res,[f(t),diff(f(t),t)]));
p:=proc(a,b,f0,omega) local fa,f1a,fb,f1b,i; global p1,p2,p3,p4;
   res(parameters=[f0,omega]);
   p1(_passed):=F(a);
   p2(_passed):=F1(a)+1;
   p3(_passed):=F(b);
   p4(_passed):=F1(b)-1;
   if type(procname,indexed) then cat(op(0,procname),op(procname))(_passed) else [p1,p2,p3,p4](_passed) end if
end proc;
assign(seq(cat(p,i)=subs(ii=i,(proc(a,b,f0,omega) p[ii](_passed) end proc)),i=1..4));
eval(p3); #Just checking the assignment
p(Pi/2,-Pi/2,1,1); #Checking p
p1(Pi/2,-Pi/2,1,1); #Checking p1
#Now solving f(a)=0, f'(a)=-1, f(b)=0, f'(b)=1 for a,b,f0,omega:
fsolve([p1,p2,p3,p4],[1.5,-1.5,1,1]);
res(parameters);
plots:-odeplot(res,[[t,f(t)],[t,diff(f(t),t)]],-3..3); #One solution
fsolve([p1,p2,p3,p4],[3,-3,.5,1]);
plots:-odeplot(res,[[t,f(t)],[t,diff(f(t),t)]],-3..3); #Another solution

It seems that there is a diff versus Diff problem. The number option is only used when the system is given by a procedure.
Try this instead:

DEplot3d(value(sys), [u1(t), u2(t), u3(t)], t = 0 ..0.4,[[u1(0) = 1, u2(0) = 1, u3(0) = 3]]);
# and
DEplot(value(sys), [u1(t), u2(t), u3(t)], t = 0 ..0.4,[[u1(0) = 1, u2(0) = 1, u3(0) = 3]],scene=[u1(t),u2(t)]);