In the definition of C3 you need to change 'm' to 'n':

C3 := n-> C1*sin(n*p*s)+C2*cos(n*p*s);

Also replace 'Sum' by 'add'.

Thirdly, B2x is a 5x5 matrix, so you can plot each element. To plot B2x[2,4] do
plot(B2x[2, 4], x = 0 .. Pi/2);

V:=16*Int(Int(r*sqrt(R^2-r^2),r=0..a/cos(theta)),theta=0..Pi/4);
#I didn't wait long enough for a result from this attempt:
#value(V) assuming R>a*sqrt(2),a>0;
#Now doing the inner integral first:
V1:=16*Int(r*sqrt(R^2-r^2),r=0..a/cos(theta));
value(V1) assuming R>a/cos(theta),a>0,cos(theta)>0;
expand(%);
Int(%,theta=0..Pi/4);
value(%) assuming R>a*sqrt(2),a>0;

#This was done in Maple 12

Suppose you have a procedure like

q:=proc(x) local a,b; a:=b; b:=34; a end proc;

Then try
q();
%;

This reflects the fact that in procedures you don't have full evaluation which is the normal behavior outside.
To avoid seeing 'b' do like this
q:=proc(x) local a,b; a:=b; b:=34; eval(a) end proc;

So in your example beginning with

eqti := (x,y) -> if x = y then 1 else 0 end if;
for i from 1 to 3 do: p[i] := ListToolsDotProduct([`$`(a[i,j],j=1..3)],[1,x,y]); end do;
for i from 1 to 3 do:
  for j from 1 to 3 do:
    eq[i,j] := subs
etc. etc.

you can in the procedure version just end with ...  eval(p[1]) end proc;  instead of just p[1] end proc;

Here is a straightforward approach first illustrated outside a procedure.
N:=8:
M:=Matrix(N,(i,j)->if i=j then e elif abs(i-j)=1 then h else 0 end if):
#Now redefine the exceptional elements:
M[1,2]:=h1: M[2,1]:=h1: M[N-1,N]:=h2: M[N,N-1]:=h2: M[N,N]:=e-b:
M;
#This could be put into a procedure like this:
p:=proc(N::posint,b,e,h,h1,h2) local M;
   M:=Matrix(N,(i,j)->if i=j then e elif abs(i-j)=1 then h else 0 end if);
   M[1,2]:=h1;
   M[2,1]:=h1;
   M[N-1,N]:=h2;
   M[N,N-1]:=h2;
   M[N,N]:=e-b;
   M
end proc;
p(7,b,e,h,h1,h2);
p(11,b,e,h,h1,h2);
#To have matrices with row or column size larger than 10 printed set 'rtablesize':
interface(rtablesize=20);
p(11,b,e,h,h1,h2);

Should anybody be interested at this late date there is a workaround at

http://www.mapleprimes.com/posts/134942-Patching-Maple?submit=135014#comment135014

The following 4 versions of a simple procedure (although a needlessly complicated one) all return the square of the input. The globals are not affected.

p:=proc(x) local a; assign(a=x); a^2 end proc;
p(8);
a;
p:=proc(x) local a; a:=x; a^2 end proc;
p(8);
a;
p:=proc(x) local a,b; a:=b; assign(a,x); b^2 end proc;
p(8);
a,b;
p:=proc(x) local a,b; a:=b; a:=x; a^2 end proc;
p(8);
a,b;

#But maybe I didn't understand your problem?

Have you looked at NonlinearFit from the Statistics package?

You are missing a semicolon after x:=cos(theta). So the space is unfortunately interpreted as an implied multiplication sign. This couldn't happen if you used one-dimensional input (aka 'Maple Input'). Compare the following:

restart;
x:= cos(Θ)*with(orthopoly); P(5,x);  
restart;
x:= cos(Θ); with(orthopoly); P(5,x);

I don't think it is a bug. The same behavior is found in Maple 12 and Maple 15 (and presumably in 13 and 14).

It is disastrous if   -1.03^b is confused with (-1.03)^b. Admittedly if the user inputs -1.03^b or -(1.03)^b then the printing is different, but you still have a product just as in your example x2:

-(1.03^b);
whattype(%);
-1.03^b;
whattype(%);

x2 := -(a)^(b):
eval(x2,[a=103/100]):
evalf(%);

whattype(%);
x3 := (-a)^(b):
eval(x3,[a=103/100]):
evalf(%);

x4 := -(a)^Pi:
eval(x4,[a=103/100]):
evalf(%);

x5 := (-a)^Pi:
eval(x5,[a=103/100]):
evalf(%);

x6 := -(a)^b:
eval(x6,[a=-103/100]):
evalf(%);

#Also try this:
evalf(-a);
#This makes the following evaluate to a float:
eval(%,a=5/6);
#Just as this does:
evalf(2*a);
eval(%,a=5/6);
#But evalf(a) just returns a. It has to. It is not a product of a type numeric and a name as -a and 2*a both are.

Try

A:=simplify( your integral )  assuming m::integer;

And after that:

simplify(A) assuming m::even;

simplify(A) assuming m::odd;

From what you call 'the required answer' it seems that your 'n' in K should have been 2*n.

For the ode approach I have a couple of suggestions. One is to keep the original ode, i.e. don't write it as a real system. Maple can handle it. The other is to use odeplot, which has the advantage of being very fast, I suppose partly because it knows to stop when hitting a singularity (and then issues a warning), whereas plot just ignores errors from dsolve and goes on to the next point which also results in an error etc.

Following Robert Lopez, let P be your polynomial. Then the following gives you pictures and animations very fast.
I use fsolve as in your worksheet to find initial conditions, but you could change that to Robert Lopez' RootOf version:

p := proc(x) option remember; global P,v,l2; local res;
if not type(x,numeric) then return 'procname'(_passed) end if;
res:=fsolve(eval(P,v=x),l2,complex);
if type(procname,indexed) then res[op(procname)] else res end if
end proc:

f:=unapply(P,l2,v):
ode:=diff(l2(v), v) = subs(l2 = l2(v), -(diff(f(l2, v), v))/(diff(f(l2, v), l2))):
sol:=(w,a)->dsolve({ode,l2(w)=p[a](w)},numeric,range=3405..5054);
#Checking that sol works:
sol(4000,1)(3500);
#Checking odeplot:
plots:-odeplot(sol(4000,1),[Re(l2(v)),Im(l2(v))],3405..5054);
#Your colors:
cols:=[red,blue,black,green,"Violet",brown,cyan,"DarkOrchid"]:
#Initial conditions given by p(4000):
plts:=seq(plots:-odeplot(sol(4000,k),[Re(l2(v)),Im(l2(v))],3405..5054,color=cols[k],thickness=3,refine=1),k=1..8):
plots:-display(plts);
#Initial conditions given by p(3500):
plts:=seq(plots:-odeplot(sol(3500,k),[Re(l2(v)),Im(l2(v))],3405..5054,color=cols[k],thickness=3,refine=1),k=1..8):
plots:-display(plts);
#Initial conditions given by p(5000):
plts:=seq(plots:-odeplot(sol(5000,k),[Re(l2(v)),Im(l2(v))],3405..5054,color=cols[k],thickness=3,refine=1),k=1..8):
plots:-display(plts);
#Animations
plots:-odeplot(sol(4000,2),[Re(l2(v)),Im(l2(v))],3405..5054,color=cols[2],thickness=3,frames=100);
plts:=seq(plots:-odeplot(sol(4000,k),[Re(l2(v)),Im(l2(v))],3405..5054,view=[-0.45..0.65,-1.3..1.3],color=cols[k],thickness=3,numpoints=500,frames=100),k=1..8):
plots:-display(plts);
plts:=seq(plots:-odeplot(sol(3500,k),[Re(l2(v)),Im(l2(v))],3405..5054,view=[-0.45..0.65,-1.3..1.3],color=cols[k],thickness=3,numpoints=500,frames=100),k=1..8):
plots:-display(plts);

#######################

Comment on the dependence on initial conditions.

Examining the simple equation x^2 + y^2 = 1 in the very same way as done in your much more complicated case exhibits clearly why parts of solutions are missing when using the ode approach: Singularities are met with.
The following mimics as much as possible the solution above.

restart;
eq:=x^2+y^2=1;
ode:=diff(subs(y=y(x),eq),x);
dsolve({ode,y(0)=1});
p := proc(xx) option remember; global eq,x,y; local res;
if not type(xx,numeric) then return 'procname'(_passed) end if;
res:=fsolve(eval(eq,x=xx),y,complex);
if type(procname,indexed) then res[op(procname)] else res end if
end proc:
p(2);
p[1](2);
plots:-complexplot([seq(p[i](x),i=1..2)], x = -3..3, color =[red,blue],thickness=3);
plots:-complexplot([-sqrt(1-x^2),sqrt(1-x^2)],x=0..3,thickness=3,color=[red,blue]);
sol:=(x,i)->dsolve({ode,y(x)=p[i](x)},numeric,range=-3..3);
sol(2,1)(.3);
plots:-odeplot(sol(2,1),[Re(y(x)),Im(y(x))],-3..3,thickness=3,refine=1);
cols:=[red,blue]:
plts:=seq(plots:-odeplot(sol(2,k),[Re(y(x)),Im(y(x))],-3..3,color=cols[k],thickness=3,refine=1),k=1..2):
plots:-display(plts);
plts:=seq(plots:-odeplot(sol(.9,k),[Re(y(x)),Im(y(x))],-3..3,color=cols[k],thickness=3,refine=1),k=1..2):
plots:-display(plts);
plts:=seq(plots:-odeplot(sol(0,k),[Re(y(x)),Im(y(x))],-3..3,view=[-1..1,-1..1],color=cols[k],thickness=3,numpoints=500,frames=100),k=1..2):
plots:-display(plts);
plts:=seq(plots:-odeplot(sol(2,k),[Re(y(x)),Im(y(x))],-3..3,view=[-1..1,-3..3],color=cols[k],thickness=3,numpoints=500,frames=100),k=1..2):
plots:-display(plts);
plots:-odeplot(sol(2,1),[x,Re(y(x)),Im(y(x))],-3..3,axes=normal,frames=100,view=[-3..3,-1..1,-3..0.5]);
plots:-odeplot(sol(2,1),[x,Re(y(x)),Im(y(x))],1..3,axes=normal,frames=100,view=[-3..3,-1..1,-3..0.5]);

Change 'pi' to 'Pi'. Then it works.

As observed by Markiyan Hirnyk in a related question the sum a+b satisfies the heat equation but with time reversed, thus has the properties of solutions of that equation when time is run backwards.
So consider that PDE (a+b will be denoted by a):

PDE:=eval(PDE1+PDE2,b=0);
#with initial and boundary conditions:
ICa:=remove(has,IC,b);
SOLa:=pdsolve(PDE,ICa,numeric,time=t,range=0..1,spacestep=1/200);
#Now suppose a solution existed on some time interval [0, t0] where t0>0 and x is in [0,1] .
Then by running the time backwards from t0 the (unique!) solution to our PDE would exist and be smooth for t in (-infinity, t0) and x in (0,1). But it is not even continuous at t = 0.
The problems will also show up numerically even if the boundary conditions are changed to allow continuity for t=0 on the x-interval [0,1].

Going backwards in time is no problem:
SOLa:-plot3d(x=0..1,t=-3..0,axes=normal);

#An additional comment: You had a timestep of 50 (!), yet consider t = 0..1. Did you mean timestep = 0.005?

restart;
with(LinearAlgebra):
#x:=2:
g:=Matrix([[3*ξ+5*x,8*x-7*ξ,2*ξ],[ξ,4*ξ+9*x,7*ξ],[5*ξ,ξ,4*x-6*ξ]]):
h:=RandomMatrix(3):
p:=g.h:
#You could make a function of x as follows:
E00:=unapply(map(int,p,ξ=1..3),x):
#Then use it like this:
E00(x);
E00(2);
E00(16);

I don't see any way of simplifying the solution as polynomials of degree 5 and above cannot in general be solved in terms of radicals. Your polynomial does have a special structure, so there could be hope. However, experimentation with randomly chosen parameter values (here integers between 1 and 9) suggests that there is no hope of a general formula.

restart;
test := [[lambda[h] = RootOf((n^2*tau[h]^2-2*n^2*tau[f]*tau[h]+n^2*tau[f]^2)*_Z^5+(-8*a*n^2*tau[h]+8*L*b*n^2*tau[h]-4*n*tau[f]*L*b+tau[h]^2*n+5*n^2*tau[f]*tau[h]+4*L*b*n*tau[h]+n^3*tau[h]^2-5*n^2*tau[f]^2+8*n^2*tau[f]*a-n^3*tau[f]^2-n*tau[f]^2-8*n^2*tau[f]*L*b)*_Z^4+(-32*a*n^2*b*L+16*n*b^2*L^2+20*n^2*tau[f]*L*b-2*n^3*tau[h]^2+tau[h]*L*b-3*tau[f]*n*a-2*tau[h]^2*n+16*n^2*b^2*L^2-12*L*b*n^2*tau[h]-4*n^2*tau[f]*tau[h]+12*n*tau[f]*L*b-16*n*b*L*a-4*L*b*n*tau[h]+tau[f]*L*b-3*a*n*tau[h]+16*a^2*n^2+7*n^2*tau[f]^2+2*n^3*tau[f]^2+10*a*n^2*tau[h]+4*b^2*L^2-3*n^2*tau[h]^2-22*n^2*tau[f]*a+2*n*tau[f]^2)*_Z^3+(-6*b^2*L^2-24*n^2*b^2*L^2-24*n*b^2*L^2-n*tau[f]^2+n^3*tau[h]^2+24*n*b*L*a+5*a*n*tau[h]-16*n^2*tau[f]*L*b-n^3*tau[f]^2-12*n*tau[f]*L*b-3*n^2*tau[f]^2-24*a^2*n^2+48*a*n^2*b*L+4*tau[f]*n*a+2*n^2*tau[h]^2+n^2*tau[f]*tau[h]-tau[h]*L*b+4*L*b*n^2*tau[h]+tau[h]^2*n+18*n^2*tau[f]*a-2*tau[f]*L*b)*_Z^2+(-tau[f]*n*a-8*n*b*L*a-16*a*n^2*b*L+tau[f]*L*b-2*a^2*n+8*a^2*n^2-2*a*n^2*tau[h]+8*n^2*b^2*L^2-2*a*n*tau[h]-4*n^2*tau[f]*a+2*b^2*L^2+8*n*b^2*L^2+4*n^2*tau[f]*L*b+4*n*tau[f]*L*b)*_Z+a^2*n)]]:
R:=rhs(test[1,1]);
p:=op(R): #The polynomial
#There is at least one real root between 0 and 1:
eval(p,_Z=0);
eval(p,_Z=1);
nms:=[op(indets(R,name))]; #The list of parameters
param1:= nms=~[1,2,3,4,5,6]; #Initial experiment
R1:=eval(R,param1);
allvalues(R1);
evalf(%);
r:=rand(1..9): #Chooses random integers between 1 and 9
param:= nms=~[seq(r(),k=1..6)];
R1:=eval(R,param);
allvalues(R1);
evalf(%);
#This loop shows that mostly no solution in terms of radicals are found:
to 10 do
   param:= nms=~[seq(r(),k=1..6)];
   R1:=eval(R,param);
   res:=allvalues(R1);
   if has([res],RootOf) then res:=evalf(res) else res end if;
   print(res)
end do:
#Giving up finding exact roots and using fsolve to find roots between 0 and 1:
to 100 do
   param:= nms=~[seq(r(),k=1..6)];
   F:=eval(op(R),param);
   res:=fsolve(F,_Z=0..1);
   if res=NULL then print("None") else print(res) end if
end do: