Here is a way using the parameters option and output=listprocedure.

restart;
f1:=(k,t)->sin(k*t);
f2:=(k,t)->cos(k*t);
g1:=(k,t)->sin(sqrt(k)*t);
g2:=(k,t)->sin(t)/(1+k+t);
sys:={diff(x(t),t)=f1(k,t)*x(t)+f2(k,t)*y(t), diff(y(t),t)=g1(k,t)*x(t)+g2(k,t)*y(t),x(0)=1,y(0)=1};
res:=dsolve(sys,numeric,parameters=[k],output=listprocedure);
X,Y:=op(subs(res,[x(t),y(t)]));
Q:=proc(k) if not k::numeric then 'procname(k)' else res(parameters=[k]); X end if end proc;
plot(Q(8)(t),t=0..2);
plot3d(Q(k)(t),k=0..10,t=0..2);
##Actually the following version is slightly faster:
Q:=proc(k,t) if not k::numeric then 'procname(k,t)' else res(parameters=[k]); X(t) end if end proc;
plot3d(Q(k,t),k=0..10,t=0..2);
#However, adding option remember to the first version makes that considerably faster:
Q:=proc(k) option remember; if not k::numeric then 'procname(k)' else res(parameters=[k]); X end if end proc;
plot3d(Q(k)(t),k=0..10,t=0..2);

You are not likely to get any kind of formula.
Consider your setup, where I corrected a couple of syntactical errors: square brackets replaced by () and comma replaced by period:

E:=5.37*theta12(t)=C*(diff(theta12(t),[t$2])+diff(theta10(t),[t$2]))+247.2*(diff(theta10(t),[t$2])*(-6.53*cos(theta12(t))+8.51*sin(theta12(t)))+diff(theta10(t),t)*diff(theta12(t),t)*(6.53*sin(theta12(t))+8.51*cos(theta12(t))))-247.2*diff(theta10(t),t)*(-diff(theta12(t),t)+diff(theta10(t),t))*(-6.53*sin(theta12(t))-8.51*cos(theta12(t)));
ics:=theta10(0)=0,theta12(0)=0, D(theta10)(0)=0,D(theta12)(0)=0;
#No success even without initial conditions for theta12
dsolve(E,theta12(t));
#Now try with a given function theta10 satisfying the initial condtions for that function:
E1:=eval(E,theta10(t)=t^2);
dsolve(E1); #No success
#However, numerical solution can be done with the given theta10 and with a value for C:
res:=dsolve({eval(E1,C=10.123456789),theta12(0)=0, D(theta12)(0)=0},numeric);
plots:-odeplot(res,[t,theta12(t)],0..2);

What is the value of k_1? Is it to be taken from your list L?
In that case by taking the first value in L you can do
res:=dsolve({eval(eq1,k_1=0.01),bcs1},numeric,method=bvp[midrich],maxmesh=1024);
plots:-odeplot(res,[eta,f(eta)],0..N);
I'm assuming that Maple uses some kind of finite difference method.

For solving the system:

res:=dsolve(eval({eq1,eq3,bcs1,bcs3},k_1=0.01),numeric,method=bvp[midrich],maxmesh=1024);
plots:-odeplot(res,[eta,theta(eta)],0..N);

Don't use the deprecated 'array', use 'Array' (capital A).

I realize that you are not using heat1 for any computations. Were you to use it for (almost) any purpose you should definitely not define it as a function.
Try with your definition
heat1(5,7); #Error
heat1(y,tau); #OK as long as y and tau evaluate to names

I found no problem with

plot(Im(x^x),x=-3..0);
plot(Re(x^x),x=-3..0);
plots:-complexplot(x^x,x=-3..0);
plot(-abs(x^x),x=-3..0);

For numerical solution infinity needs to be finite!
In this example infinity is N = 20.
Also the singularity at x=0 is avoided by replacing 0 by eps.
Thus the result is highly questionable.
But you can play with it.

restart;
ode:=diff(y(x),x,x)+diff(y(x),x)/x-(1-1/x^2)*sin(2*y(x))/2+0.3*sin(y(x))=0;
#This direct approach doesn't work:
res:=dsolve({ode,y(0)=Pi,y(10)=0},numeric,method=bvp[midrich]);
#Shooting from x = N = 20:
N:=20: #"Infinity"
respar:=dsolve({ode,D(y)(N)=y1,y(N)=0},numeric,parameters=[y1],output=listprocedure);
Y:=subs(respar,y(x));
eps:=1e-7: #"Zero"
Q:=proc(y1) respar(parameters=[y1]); Y(eps)-evalf(Pi) end proc;
Q(.1);
Q(-.1);
Q(0);
fsolve(Q,0);
respar(parameters=[%]);
plots:-odeplot(respar,[x,y(x)],eps..N);

Remark.
Be aware that there are (at least) 2 solutions to Q=0, one negative, the other positive.
plot(Q,-.05..0.05,adaptive=false);
#So make fsolve pick the positive one
fsolve(Q,0.002);

I picked the constants a,b,...,h arbitrarily:
eq1:=a*diff(phi(x),x$4)+b*diff(phi(x),x,x)+c*phi(x)+d*diff(psi(x),x,x)+e*psi(x)=0;
eq2:=d*diff(phi(x),x,x)+e*phi(x)+j*diff(psi(x),x,x)+h*psi(x)=0;
bcs:=phi(0.83)=1,D(phi)(0.83)=0,phi(-0.83)=1,D(phi)(-0.83)=0,psi(0.83)=1,psi(-0.83)=1;
res:=dsolve(eval({eq1,eq2,bcs},{a=.1,b=.2,c=.3,d=.1,e=.2,j=.3,h=.4}),numeric);
plots:-odeplot(res,[[x,phi(x)],[x,psi(x)]],-0.83..0.83);

Maple has manuals. In Maple go to Help/Manuals, Resources, and more/Manuals/Programming Guide.

A pdf version can be found at
http://www.maplesoft.com/documentation_center/

Assuming that R is a given function of P you can do as follows:

sys:=diff(y(x),x) +y(x)^2 =P(x), diff(P(x),x) = R(P(x));
#First example
R:=sin:
res1:=dsolve({sys,y(0)=0,P(0)=1},numeric);
plots:-odeplot(res1,[[x,y(x)],[x,P(x)]],-15..10);
#Second example
R:=p->p/(1+p^2):
res2:=dsolve({sys,y(0)=0,P(0)=1},numeric);
plots:-odeplot(res2,[[x,y(x)],[x,P(x)]],-15..10);

You don't provide an explicit context, so I took an example out of my head:

restart;
ode1:=diff(x(t),t)=sin(x(t)-sin(t));
res1:=dsolve({ode1,x(0)=0},numeric,output=listprocedure);
X1:=subs(res1,x(t));
plot(X1,0..50); #You can also use plot(X1(t),t=0..50);
ode2:=diff(x(t),t)=cos(x(t)-sin(t));
res2:=dsolve({ode2,x(0)=0},numeric,output=listprocedure);
X2:=subs(res2,x(t));
plot(X2,0..50);
plot(X2-X1,0..50);
plot([X1,X2,sin],0..50);
plot([X1(t),X2(t),sin(t)],t=0..50);
##############
If you have a system you can do differently. Here (artificially) we make our two euations into one system:
res:=dsolve({ode1,subs(x=y,ode2),x(0)=0,y(0)=0},numeric);
plots:-odeplot(res,[[t,y(t)-x(t)],[t,sin(t)]],0..50,caption="Difference y(t)-x(t) and sin(t)");

If you try L := [0., .2, .43]; you will see that your code runs. Problems start at A=0.44.

You could look into continuation in A. See
?dsolve,numeric,bvp

Continuation here would mean replacing A by A0*(1-c)+A1*c in the equations, where A0 could be 0.2 and A1 could be 0.5. Then use 'continuation=c' in dsolve. It is definitely not guaranteed to work, but you could experiment with it.

In the help page for collect it says about the optional argument func:
"A function may be specified using the optional argument func. It is applied to the coefficients of the collected result. Often simplify or factor will be used."

So in the case you encountered the simplifier used was normal.

@snowww You need to be aware that the optional argument to pdsolve, 'build', does not in general give you the general solution. It just gives you building blocks.
Thus you cannot expect to make a single 'building block' fit your boundary conditions.
In general you will need an infinite series of such building blocks, just as is the case for the usual textbook examples for the heat equation.

pde:=(1-x^2)*diff(c(x,y),y)=diff(c(x,y),x,x)-alpha^2*c(x,y);
res:=pdsolve(pde); #Finding solutions by separating variables
pdsolve(pde,build); #solving the odes too
op(2,res);
dsolve(op([2,1,1],res));
dsolve(op([2,1,2],res));

#You could use numerical solution instead. You need to make alpha concrete. I have taken alpha=1:

ibcs:={c(x,0)=1,c(0,y)=1,D[1](c)(1,y)=0};
sol:=pdsolve(eval(pde,alpha=1),ibcs,numeric);
sol:-plot3d(x=0..1,y=0..1);

To answer the question: The locals are names.
From the help page for 'procedure':

"Local variables that appear in the local localSequence; clause may optionally be followed by :: and a type. As in the case of the optional returnType, this is not a type declaration, but rather an assertion. If kernelopts(assertlevel) is set to 2, any assignment to a variable with a type assertion is checked before the assignment is carried out. If the assignment violates the assertion, then an exception is raised."

Try setting assertlevel:
kernelopts(assertlevel=2);
Then run your procedure
Max2(1,7,5,4,6,10,35,63.5,-10,5);
Error, (in Max2) assertion failed in assignment, expected string, got 1
#Another example:
p:=proc(x) local a::float; a:=x; a end proc;
p(6);
p(7.);

A simple Google search for "bubble sort maple" brought up

http://www.mhhe.com/math/advmath/rosen/r5/student/ch09/maple.html

by Kenneth H. Rosen, AT&T Laboratories (with very minor changes by me: array->Array).

BubbleSort:=proc(L::list)
  local i, j, temp, T;
  T:= Array(L);
  for i from 1 to nops(L)-1 do
    for j from 1 to nops(L)-i do
      if T[j] > T[j+1] then
        temp:=T[j];
        T[j] := T[j+1];
        T[j+1] := temp;
      end if
    end do
  end do;
  convert(T,list)
end proc:

BubbleSort([3,2,4,1,5]);