The first thing to notice is that j and k are determined uniquely by the requiremens n=2^j+k and 0<=k<=2^j-1.
In fact then j = ilog2(n) (a Maple function).
Thus you can do:
restart;
H:=t->piecewise(t<0,0,t<1/2,1,t<1,-1,0);
h:=proc(n::nonnegint,t) local j,k; if n=0 then return piecewise(t<0,0,t<1,1,0) end if;
  j:=ilog2(n);
  k:=n-2^j;
  2^(j/2)*H(2^j*t-k)
end proc;
h(0,t);
h(1,t);
S:=[seq(plot(h(n,t),t=0..1,caption=typeset("n= ", n)),n=0..31)];
plots:-display(S,insequence);

Let me present the easiest solution first.
restart;
ode1 := (K[Q]*T*R[a]/K[kol]+R[a]*B[m]*sqrt(T/K[kol]))/K[i]+K[b]*sqrt(T/K[kol]) = 0;
solve(ode1,T);
#There are two solutions for T, we choose the nontrivial one:
ode1:=T=%[2];
ode2 := (1/2)*(-(4*(diff(theta(t), t)+theta(t)))*M+2*B1*(diff(theta(t), t)+theta(t))-2*w[2]*sin(theta(t))+m1*g*sin(theta(t))*l[kol]-2*B1*theta(t)+2*w[1]*sin(theta(t))-2*m1*g*sin(theta(t))*l[2])/l[1];
ode3 := subs(T = (1/2)*(-(4*(diff(theta(t), t)+theta(t)))*M+2*B1*(diff(theta(t), t)+theta(t))-2*w[2]*sin(theta(t))+m1*g*sin(theta(t))*l[kol]-2*B1*theta(t)+2*w[1]*sin(theta(t))-2*m1*g*sin(theta(t))*l[2])/l[1], ode1);
##Now you simply work with ode3 for which diff(theta(t),t) can easily be isolated by convertsys (and us).
##Assign constants as you did. Then
dsol:=dsolve({ode3,ic1},numeric);
odeplot(dsol,[t,theta(t)],0..10);
##############################
##The messier solution does essentially the same thing.
The problem is that there are two solutions for the derivative of theta when trying to isolate diff(theta(t),t) from eol as convertsys must.

You can do that yourself and then decide which is intended.
Before you assign values to the constants do the following:
eqs:=solve(eol=0,{diff(theta(t),t)});
#Notive that eqs is a sequence of 2 sets.
#Checking the right hand sides for theta(t)=0 for all t.
eval((rhs@op)~([eqs]),theta(t)=0);
#We see that eqs[1] has theta=0 as an equilibrium point. This means that if theta(0)=0 then theta(t) = t for all t.
#Thus you are probably more interested in the other solution.
#Now assign values to the constants as you do in your worksheet.
#Then consider the second solution:
dsol1b:=dsolve(eqs[2] union {ic1},numeric);
plots[odeplot](dsol1b, [t, theta(t)], 0 .. 10);
#It appears that theta(t) approaches an equilibrium solution.
#We check to see if there is an equilibrium for eqs[2] near theta=-4000:
subs(theta(t)=theta,rhs(op(eqs[2])));
fsolve(%=0,theta=-4000);
#Indeed!

I asked myself the same question several years ago and finally came up with the following procedure, which I called NestDo.

NestDo:=proc() local X1,X,N,null,R,u,k,K,ff,j,i,tmp,SEQ,expr,nvn,a,ud,kk,kmd,SUBS,sts,A;
X1:=[_passed];
N:=nops(X1);
X:=NULL;
for j from 1 to N do
   if nops(X1[j])=1 then X:=X,[op(X1[j]),null] else X:=X,X1[j] end if;
   if not member(nops(X1[j]),{1,2}) or not type(op(1,X1[j]),{name={range,list,set},range,list,set}) then error "Each list must contain 1 or 2 elements, where the first has type range, list, set or name={range,list,set}."
   end if;
end do;
X:=[X];
if foldl(`and`,true,seq(type(op(1,X[j]),name={range,list,set}),j=1..N)) then
    expr:=true
  elif foldl(`and`,true,seq(type(op(1,X[j]),{range,list,set}),j=1..N)) then
    expr:=false;
  else
    error "The lists must all have the type name={range,list, set} or all have the type range,list or set."
end if;  
R:=map2(op,1,X);
if not expr then R:=zip(`=`,[seq(cat(nvn,j),j=1..N)],R) end if;
if expr then
    k:=map(lhs,R)
  else
    k:=[seq(cat(nvn,j),j=1..N)]
end if;
K:=[seq(k[j..N],j=1..N)];
##
##The actual commands (sts) are saved and initially replaced with space holders (kmd):
##
sts:=map2(op,2,X);
if expr then
   for j from 1 to N do A[j]:=unapply(sts[j],K[j]) end do;
   sts:=[seq(A[j],j=1..N)]
end if;
kmd:=[seq(cat(a,N+1-nops(i))(op(i)),i=K)];
u:=ListTools:-Reverse(kmd);
tmp:=foldl(ff,u,seq(R[j],j=1..N));
SEQ:=subs(ff=seq,tmp);
u:=ListTools:-Flatten(eval([SEQ]));
u:=ListTools:-MakeUnique(u);
SUBS:=zip(`=`,[seq(cat(a,j),j=1..N)],sts);
ud:=NULL;
for kk in u do ud:=ud,eval(kk,{op(SUBS)}) end do;
eval([ud],null=NULL)
end proc:
##############
Notice that the inputs are lists containing information about range and action (if any) at each level.
The inner loop is the first, the outer the last.
## A few examples
res:=NestDo([k1=1..2,f(k1,k2,k3)],[k2=1..3,g(k2,k3)],[k3=1..3,h(k3)]);
#Compare with
printlevel:=3:
for k3 from 1 to 3 do
  h(k3);
  for k2 from 1 to 3 do
      g(k2,k3);
    for k1 from 1 to 2 do
      f(k1,k2,k3)
    end do
  end do
end do;
printlevel:=1:
k1:='k1': k2:='k2': k3:='k3':
NestDo([k1=1..k2,g(k1,k2)],[k2=1..3]);
#The cartesian product in Carl's example:
NestDo([k1=[e,f],[k3,k2,k1]],[k2=[c,d]],[k3=[a,b]]);

You could start by visualizing the problem. Here is an animation with k as the animation parameter:

p:= x^2+2*x+3;
q:= k+5*x-7*x^2;
plots:-animate(plot,[[p,q],x=-2..2],k=1..25);

You see that you will have to find the x-values of the points of intersection between the two curves.
Use solve for that. You will find two values r1 and r2 (say).
Now find the area A by integrating q-p from x=r1 to x=r2 (assuming that r1<r2).
The area A will depend on k.
Use solve to find k from the equation A=36.

You are obviously interested in nontrivial solutions!

I don't think there is any way you can make pdsolve solve this problem analytically.
The reason is found if you try using separation of variables as you started doing.

restart;
pde := IE*(diff(w(x, t), x, x, x, x))+m*(diff(w(x, t), t, t)) = 0;
s := pdsolve(pde, w(x, t));
ode1 := op([2, 1, 1], s);
ode2 := op([2, 1, 2], s);
##We want to make the solutions to ode2 satisfy the boundary conditions.
## We first consider the case _c[1] >=0 and make the following replacement of _c[1]:
ode2p:=eval(ode2,{_c[1]=c^4*L^(-4),_F1=F});
sol2 := dsolve(ode2p);
#Using the boundary conditions:
eqs:=[seq(eval(diff(rhs(sol2),[x$k]),x=0)=0,k=0..1),seq(eval(diff(rhs(sol2),[x$k]),x=L)=0,k=2..3)];
LinearAlgebra:-GenerateMatrix(eqs,[_C1,_C2,_C3,_C4]);
A:=%[1];
#For nontrivial solutions to exist we need the determinant of A to be zero:
dt:=LinearAlgebra:-Determinant(A);
dt1:=factor(combine(dt));
For c=0 we have the solution  
eval(sol2,c=0);
# i.e. a constant.
#Other c-solutions:
dt1:=simplify(dt1*c^(-6)*L^6);
dt2:=collect(dt1,cos,x->convert(x,cosh));
eq2:=cos(c)=-sech(c);
#The equation eq2 cannot be solved analytically so that is the problem!
solve(eq2,c);
## Plotting
plot([lhs,rhs]~(eq2),c=0..15);
##We can easily find the roots numerically since we see that they must be close to the zeros of cosine.
# Even for moderately high values of k they are very close indeed.
Digits:=15:
for k from 0 to 6 do r[k]:=fsolve(eq2,c=(2*k+1)*Pi/2) end do;
evalf(seq((2*k+1)*Pi/2,k=0..6));
####
##One should consider if it is possible that _c[1] could be negative. If you try replacing  _c[1] by -c^4*L^(-4) (with c > 0) and proceed as above you will find that the determinant is never zero. 

###########################Plotting solutions corresponding to each eigenfrequency########
##Example with c = r[1] ####
##Before evaluating the roots r[k] for k=1..6 notice that I set Digits:=15:
ode1;
eval(ode1,{_c[1]=c^4*L^(-4),_F2=G});
ode1a:=subs(IE=K*m,%);
sol1:=dsolve(ode1a);
sol2;
A1:=eval(A,{c=r[1],L=1});
LinearAlgebra:-Determinant(A1); #Just a check
LinearAlgebra:-GaussianElimination(A1);
A10:=fnormal~(%); #Roundoff removed
C:=LinearAlgebra:-LinearSolve(A10,<0,0,0,0>,free=p);
C:=eval(C,p[1]=1);
param2:={_C1=C[1],_C2=C[2],_C3=C[3],_C4=C[4],L=1,c=r[1]};
param1:={_C1=3,_C2=-2,L=1,K=1,c=r[1]};
S:=eval(rhs(sol2),param2)*eval(rhs(sol1),param1);
plot3d(S,x=0..1,t=0..1);
plots:-animate(plot,[S,x=0..1],t=0..1,frames=100);
########### The other eigenfrequencies ###
## A quick and dirty plotting procedure relying on the results above.
Psol:=proc(k::posint,T::positive:=1/k^2) local Ak,C,p,param1,param2,S;
   uses LinearAlgebra;
   Ak:=eval(A,{c=r[k],L=1});
   Ak:=fnormal~(GaussianElimination(Ak));
   C:=eval(LinearSolve(Ak,<0,0,0,0>,free=p),p[1]=1);
   param2:={_C1=C[1],_C2=C[2],_C3=C[3],_C4=C[4],L=1,c=r[k]};
   param1:={_C1=3,_C2=-2,L=1,K=1,c=r[k]};
   S:=eval(rhs(sol2),param2)*eval(rhs(sol1),param1);
   plots:-animate(plot,[S,x=0..1],t=0..T,frames=100)
end proc:
Psol(1);
Psol(2);
Psol(3);
Psol(4);

To show that x is in the column space of B you need only show that the equation B.z = x has a solution for z, since B.z is just a linear combination of the columns of B, the coefficients being z[1], z[2], ... .
The row space of A is the column space of the transpose A^%T so in your situation you can do:

A:=<<2,4,1>|<1,-1,1>|<3,3,2>>;
b:=<8,5,4>;
x:=LinearAlgebra:-MatrixInverse(A,method=pseudo).b;
A^%T.z=x; #The equation to solve by LinearSolve below
LinearAlgebra:-LinearSolve(A^%T,x,free=t);

#And yes indeed there is a solution, and if there is one there is necessarily infinitely many, since A is not invertible. Thus the free parameter.

Since you only have to determine if there is a solution (not what it is) you can just do

LinearAlgebra:-GaussianElimination(<A^%T|x>);

and observe that the ranks of the augmented matrix <A^%T|x> and A^%T are equal.

dsolve simply cannot find a formula for the solution.

You can solve numerically. Let your ode be called 'ode'.

solnum:=dsolve({ode ,y(-2)=1.3,D(y)(-2)=0.7},numeric);
plots:-odeplot(solnum,[x,y(x)],-2..1);

To understand the complexity of the solution obtained by the plain use of dsolve as Carl mentions in his answer, try

ode:=EJ*diff(w(x),x,x,x,x) - (Nr/2*(1+x)+Nl/2*(1-x))*diff(w(x),x,x) = q(x);
ode2:=subs(diff(w(x),x,x)=w2(x),ode); #Second order
res2:=dsolve(eval(ode2,q=0)); #For simplicity consider the homogeneous equation
convert(res2,Bessel); #AiryAi and AiryBi can be converted to Bessel functions

The definition of b^alpha in general is

b^alpha = exp(alpha*ln(b));

where ln(b) is the principal value of ln, i.e, the value of ln that has argument in the half open interval (-Pi,Pi].
Certainly in Maple the principal value of arg is used.

Recall that ln(b) = ln(abs(b)) + I*arg(b) where this time the ln in ln(abs(b)) is just your usual logarithm defined on the positive reals.

Thus you have:

alpha:=1/10; #On purpose I avoided floats, since that is not the issue.
b:=-1/5;
b^alpha;
evalc(%);
evalf(%);
b^alpha=exp(alpha*ln(b));
evalc(%);
ln(b);

restart;
alias( seq(z[i]=z[i](seq(z[j],j=3..4)),i=1..2));
alias(f=f(seq(z[j],j=3..4)));
alias(seq(seq(Dp[m,n]=diff(z[m](seq(z[j],j=3..4)),z[n]),m=1..2),n=3..4));
f:=z[1]-z[2]*z[3];
for j from 3 to 4 do
df_dz[j]:=diff(f,z[j]);
od;
df_dz[3];
convert(df_dz[3],D);
subs(z[2]=4711,%); #Result without removing aliases
##Now removing aliases
alias(alias()=~alias()); #Removing all existing aliases
##More explicitly like this:
##alias(seq(z[i]=z[i],i=1..2), f=f, seq(seq(Dp[m,n]=Dp[m,n],m=1..2),n=3..4));
df_dz[3];
convert(%,D);
subs(z[2](z[3],z[4])=4711,%);

############
##Comparing with version without aliases:
restart;
f:=z[1](z[3],z[4])-z[2](z[3],z[4])*z[3];
f:=unapply(f,z[3],z[4]);
for j from 3 to 4 do
df_dz[j]:=diff(f(z[3],z[4]),z[j]);
od;
df_dz[3];
convert(df_dz[3],D);
subs(z[2](z[3],z[4])=4711,%);

The help page gives the following example:

with(LinearAlgebra):
A := DiagonalMatrix([5,12,13]);
B := <<5,0,0>|<3,12,0>|<4,-4,13>>;
P:= IsSimilar(A, B, output=['C']); #I omitted 'query'
A = P^(-1) . B . P;

Thus the answer is No!

Obviously you must have written sinx instead of sin(x) in your differential equation.
If I'm mistaken then post the code for de.

To Maple sinx is just a name like any other.

Try this:
if Pi>3 then "Ok" end if;
evalb(Pi>3);
is(Pi>3);
evalb(evalf(Pi)>3);

Solution:
Change x:=Pi/3;  to x:=evalf(Pi/3);

Another thing:

By n:=+2 you probably meant n:=n+2; And that is how it should be done then.

To get started you would have to make it easier to change the stepsize in your worksheet. A common method of estimating errors involves computing the solution for stepsize h and also for stepsize h/2.

I propose to let your number 10 of steps be N. That obviously involves changing 9 and 11 into N-1 and N+1, respectively.
I have done that in your worksheet and have made a few other changes such as using 'a' for the end of the interval and NOT using assign but using subs or eval in the relevant places.

Also I made a loop over N=[5,10,20] and compared graphs.

These graphs don't prove anything, but can be reassuring. They are a start.

MaplePrimes15-02-04odebvp_discrete.mw

solve will succeed if it can find a formula for the solution. If none such exists as I believe is the case here, you could try numerical solution using fsolve. For that a and b has to be made concrete, i.e. given numerical values.
Notice that the mathematical constant Pi is spelled Pi not pi.
Example:
fsolve(eval(gl1,{a=.1,b=1}),z);