Does Not(infinity) serve the purpose:

type(Pi,Not(infinity));

or maybe rather

type(Pi,And(realcons,Not(infinity))); #complexcons?
type(a,And(realcons,Not(infinity)));

Use print.

testproc:=proc() local c;
  c:=plottools:-circle([0,0],1,color=red);
  print(plots:-display(c));
  2
end proc;

The second ode is pretty simple. Integrate it once and substitute diff(sigma(y),y)) = C -diff(T(y),y) into the first:

restart;
Eq2 := diff(T(y), y, y)+(diff(sigma(y), y))*(diff(T(y), y))+(diff(T(y), y))^2+(exp(-y)+exp(y))^2 = 0;
Eq3 := diff(sigma(y), y, y)+(diff(T(y), y, y))= 0;
bcs2:=T(h1)=0,T(h2)=1,sigma(h1)=0,sigma(h2)=1;
Seq:=int(lhs(Eq3),y)=C;
subs(diff(sigma(y),y)=C-diff(T(y),y),Eq2);
Teq:=simplify(%);
dsolve({Teq,T(h1)=0,T(h2)=1});
Tres:=collect(%,[exp(2*y),exp(-2*y),exp(-C*y),y],factor);
Sres:=dsolve({eval(Seq,Tres),sigma(h1)=0});
Ceq:=eval(rhs(Sres),y=h2)=1;
Cres:=solve(Ceq,{C});
Sres2:=subs(Cres,Sres);
Tres2:=eval(Tres,Cres);
odetest({Tres2,Sres2},{Eq2,Eq3,bcs2});

Your original problem as well as the revised one in the comments doesn't have any solution.
In my comment to the revised one I pointed out that if n is not an integer (so sin(n*Pi)<>0) then the integral A=Int(F(z)^n*h(z),z=0..2), where F(z) is real-valued and h(z)>=0, cannot be zero unless F(z)*h(z) = 0 for (almost) all z in 0..2.
Using the definitions from my first answer:
M_F:=2*l*(z-l)-z^2/2;
M_1:=piecewise(z<l, 2*l*(z-l)+(l-z)/q-z^2/2, 2*l*(z-l)-z^2/2);
M_2:=2*l*(z-l)+(2*l-z)/q-z^2/2;

We get
f:=(M_F+X_1*M_1+X_2*M_2)*M_1; #This is F(z)*h(z)
simplify(f) assuming z<1;
collect(%,z); #No chance of making that identically zero, but we continue:
coeffs(%,z);
solve({%},{_X1,X_2});

The same negative result comes out in the revised version.

I had no problem at all (in Maple 18) with

restart;
ode:=diff(Theta(t), t, t)+2*Zeta*omega[balanceue]*(diff(Theta(t), t))+omega[balanceue]^2*Theta(t) = M[p]/m[balanceue];
ICS := Theta(0) = (1/8)*Pi, (D(Theta))(0) = 0;
res1:=dsolve({ICS, ode}, Theta(t), method = laplace);
res2:=dsolve({ICS, ode}, Theta(t));

Of course I cannot solve numerically without knowing the values of the constants.

Anyway using the constants
params:={Zeta=.3,omega[balanceue]=1,m[balanceue]=1,M[p]=1};

# I did:
res10:=eval(res1,params);
res20:=eval(res2,params);
plot(rhs(res10)-rhs(res20),t=0..5); #so res10=res20
res3:=dsolve(eval({ICS, ode},params),numeric);
plots:-odeplot(res3,[t,Theta(t)-rhs(res10)],0..5); #Pretty close

Your use of the word parameter is confusing (me). Presumably you mean 3 dependent variables x,y,z. The system has one parameter: r.

The system changes stability at r = 0. The system is linear and can be written u' = A.u, where u=(x(t),y(t),z(t)).
The eigenvalues of A are easily found in Maple (or by hand).
The solution itself can be found by
dsolve({sys_ode,x(0)=x0,y(0)=y0,z(0)=z0});

Another idea is to use that the sum of x, y, and z is constant, so that the 3D-system can be reduced to a 2D-system, still linear but now inhomogeneous:

sys_ode[1]+sys_ode[2]+sys_ode[3];
sys2:=subs(z(t)=x0+y0+z0-x(t)-y(t),[sys_ode[1..2]]);

You could try

res:=Optimization:-LSSolve([E[1],E[2],E[3]],C[1]=-1..1,C[2]=-1..1,C[3]=-1..1,iterationlimit=10^5);

I changed a few things. Besides syntactical changes I divided out a q factor from each equation.

restart;
B:=1:
q:=1*10^3:
l:=1:
n:=4.7:
M_F:=2*l*(z-l)-z^2/2; #No need to make the Ms functions. Divided by q
M_1:=piecewise(z<l, 2*l*(z-l)+(l-z)/q-z^2/2, 2*l*(z-l)-z^2/2); #Divided by q
M_2:=2*l*(z-l)+(2*l-z)/q-z^2/2; #Divided by q
#one_ and two_int should be functions of X_1 and X_2, not z.
one_int:=proc(X_1,X_2) local res;
  res:=evalf(Int(B*(M_F+X_1*M_1+X_2*M_2)^n*M_1,z=0..2*l));
  print([_passed],res);
  res
end proc;
two_int:=proc(X_1,X_2) local res;
  res:=evalf(Int(B*(M_F+X_1*M_1+X_2*M_2)^n*M_2,z=0..2*l));
  print([_passed],res);
  res
end proc;
#Printing in these procedures lets us know what is going on.
##The following prints out for quite a while and then stops printing, but it keeps evaluating. I interrupted the computation.
fsolve([one_int,two_int],[-20,19]);
#The last printout I got was for input:
res:=[-1.0264919886644220, 0.25832794180969000000000e-1];
one_int(op(res));
two_int(op(res));

To plot the curve given by a parametrization as here do (for the second as an example):

plot([cos(t),sin(t)*cos(t),t=0..2*Pi]);
#More fun to animate:
plots:-animate(plot,[ [cos(t),sin(t)*cos(t),t=0..T] ],T=0..2*Pi);

You cannot change the identity matrix if defined as such. But you can make A a matrix with the same elements.
However, see the comment by Kitonum and my response below.
Also certainly the final value in a for loop must be of type numeric, i.e. n has to have a concrete value.
(When your error message mentions also 'character' it refers to strings: see footnote).

restart;
n:=5: #Example
for i to n do B[i]:=LinearAlgebra:-RandomMatrix(n) end do; #Example B's.
A:=Matrix(n,(i,j)->if i=j then 1 else 0 end if);
for i from 1 to n do A:=A.B[i] end do:
A;

Alternatively use the inplace option:

A:=Matrix(n,(i,j)->if i=j then 1 else 0 end if);
for i from 1 to n do LinearAlgebra:-Multiply(A,B[i],inplace) end do:
A;

##########################
Footnote:
for k from "d" to "m" do k end do;

Rather than use the answer I would suggest using the original ode, which we can easily recover:

restart;
y(t) = _C1*exp(-1.*t)*sin(.57736*t)+_C2*exp(-1.*t)*cos(.57736*t);
#The roots of the characteristic polynomial are:
r1:=-1+0.57736*I;
r2:=conjugate(r1);
#Thus the polynomial is:
p:=simplify(expand((lambda-r1)*(lambda-r2)));
with(DEtools):
#Although it is easy to write down the ode from p here I shall use diffop2de:
ode:=diffop2de(p,y(t),[lambda,t])=0;
#A plot of 7 solution curves:
DEplot(ode,y(t),t=0..5,[seq([y(0)=y0,D(y)(0)=0],y0=-3..3)]);
#An animation, where the initial slope is the animation parameter:
plots:-animate(DEplot,[ode,y(t),t=0..5,[seq([y(0)=y0,D(y)(0)=y1],y0=-3..3)]],y1=-2..2);

I had success with an initial elimination:

restart;
eqns := [A = (gr_c+delta)*kh^(1-alpha)/sav_rate, theta = Rk*(Rh-Rk)/(gamma*((Rh-Rk)^2+sigma^2)), theta = 1*1+kh, Rk = 1+rk-delta, Rh = 1+rh-delta, rk = A*alpha*kh^(alpha-1), rh = A*(1-alpha)*kh^alpha, sigma = sigmay/theta, varrho = Rap^((eps-1)*eps/(1-gamma)), Rap = Rk^(1-gamma)+(1-gamma)*Rk^(-gamma)*theta*(Rh-Rk)-.5*Rk^(-gamma-1)*gamma*(1-gamma)*theta^2*((Rh-Rk)^2+sigma^2), R = Rk+theta*(Rh-Rk), beta = ((1+gr_c)/R)^(1/eps)/varrho];
###
vals := [alpha = .36, delta = 0.6e-1, sigmay = sqrt(0.313e-1), gamma = 3, eps = .5, gr_c = 0.2e-1, sav_rate = .23];
###
eqns_a:=eval(eqns, vals);
eliminate(eqns_a,{A,theta,Rk,Rh,sigma,varrho});
S1,S2:=op(%);
vars2:=indets(S2,name);
res1:=fsolve(S2, vars2=~(0..infinity));
res2:=eval(S1,res1);
#Check:
eval((lhs-rhs)~(eqns_a),res1 union res2);

Using Maple all the way:

restart;
pde:=diff(u(x,t),x,x)=u(x,t)+diff(u(x,t),t);
bcs:=u(0,t)=0,u(Pi,t)=0;
ic:=u(x,0)=Pi-x/2;
res:=pdsolve({pde,bcs,ic});
##You may want to replace the index _Z1 with something more familiar:
res1:=subs(_Z1=n,res);
##Whether or not you do that you can continue with an animation with a truncated series TS, here the first 20 terms:

TS:=value(eval(rhs(res1),infinity=20));
plots:-animate(plot,[TS,x=0..Pi],t=0..1);

##Actually you can use the whole series, but for t=0 you don't get a picture:
plots:-animate(plot,[rhs(res),x=0..Pi],t=0..1);
#The absence of a graph for t = 0 is clearly due to the fact that evalf cannot get a satisfactory result as in this example with x = 1.234567:
evalf(eval(rhs(res),{t=0,x=1.234567}));
#Returns unevaluated.
## You can put the initial graph in as a background:
p1:=plot(Pi-x/2,x=0..Pi);
plots:-animate(plot,[rhs(res),x=0..Pi],t=0..1,background=p1);

As I see it, it basically comes down to what you expect as output of:

int(Dirac(a)/a,a=0..100); #Maple: undefined
evalf(Int(Dirac(a)/a,a=0..100)); # Maple: 0.

Comment. You may want to try these simple examples:
A:=Int(Dirac(x-1),x=0..2);
value(A); #OK
evalf(A); #OK
evalf(Int(Dirac(x-1),x=0..2,method=_CCquad)); #OK
evalf(Int(Dirac(x-1),x=0..2,method=_NCrule)); #Returns unevaluated
evalf(Int(Dirac(x-1),x=0..2,method=_d01ajc)); #NAG routine returns Float(undefined)
#integration done by dsolve:
dsolve({diff(y(x),x)=Dirac(x-1),y(0)=0}); #OK
B:=dsolve({diff(y(x),x)=Dirac(x-1),y(0)=0},numeric);
B(2); #returns 0.

All solutions to 'ode' satisfy y(0)=0:
dsolve(ode);
eval(%,x=0);
dsolve(ode,y(x),series);
#thus the following results in NULL, i.e. no solution:
dsolve({ode,y(0)=1});

#####An illustration:
DEtools[DEplot](ode,y(x),x=0..1,[seq([y(1)=0,D(y)(1)=i/10],i=-10..10)],linecolor=gray);
DEtools[DEplot](ode,y(x),x=0..1,[seq([y(1)=-i/10,D(y)(1)=i/10],i=-10..10)],linecolor=gray);
#Notice that in both cases no matter what your conditions are at x=1 the curves will end up at x=0.
#At the risk of belaboring the point here is yet another illustration:
res:=dsolve({ode,y(1)=y0,D(y)(1)=y1},numeric,parameters=[y0,y1]);
p:=proc(y0,y1) res(parameters=[y0,y1]); plots:-odeplot(res,[x,y(x)],0..1) end proc:
plots:-animate(p,[1,y1],y1=-10..20,trace=24);