Maybe 'indets' can do what you want.

If you provided an example more detail could be given.

Again answering your latest question:

In order to make it possible to follow this I bring the whole code:

restart;
with(plots):
eq := z-> (H^2+(-1)*.27*(1+z)^3-(1/20000)*(1+z)^4)/(H^2)^.1 = .7299500000;
Y := z->if not type(z,numeric) then 'procname(z)' else fsolve(eq(z), H=1) end if:
p := dsolve({D(L)(z) = L(z)/(1+z)+(1+z)/Y(z),L(0)=0}, type=numeric, range=0..5,known=Y);
odeplot(p,[z,L(z)],0..5);

#The integral you are asking about:

R := evalf(Int(1/Y(z), z=0..10));

#If by undefined integral you meant to say indefinite integral (i.e. in this case something like G(z1)=int(1/H(z),z=0..z1) considered as a function of z1) then it is best to consider the systems approach used earlier and this time introduce a differential equation for G(z).

yp := implicitdiff(eq(z), H, z);
ode:=diff(H(z),z)=subs(H=H(z),yp);
ode2:=diff(G(z),z)=1/H(z);
p := dsolve({D(L)(z) = L(z)/(1+z)+(1+z)/H(z), ode,ode2,L(0)=0,H(0)=1,G(0)=0}, numeric, output=listprocedure);

odeplot(p,[z,G(z)],0..10);
Gp,Hp,Lp:=op(subs(p,[G(z),H(z),L(z)]));
Gp(10);

You are using G as a function, thus it must be defined as one:

G := Heaviside: #No x

int(G(x)*G(t-x), x = 0 .. 3);

If I'm correct in assuming that A has entries that are homogeneous polynomials then the problem is equivalent to asking the following question:

Given n  nxn matrices A_i.

Can we determine vectors b_i such that the sequence of equations A_i . z = b_i (i = 1..n) have a common solution vector z.

Solution: Pick any z-vector, and let b_i be determined by the equations.

As Axel Vogt suggests all int's have been replaced by Int's, and no evalf's are present. I let plot do that itself.

Re and Im have been replaced by evalc(Re(...)) and evalc(Im(..)). evalc makes the assumption that any variable name represents a real number as  in your case where t and epsilon are real. Thus assume(epsilon,real) has been removed.

In a couple of places superfluous Re's have been removed.

The plot took less than half a minute.

I have uploaded my revised version MaplePrimes10-11-0.mw

F := X-> Vector([f(X[1],X[2],X[3]),g(X[1],X[2],X[3]),h(X[1],X[2],X[3])]):

mtaylor~(F(X),[[X[1],X[2],X[3]]$3],3);
G := X-> Vector([exp(X[1])+X[2],(X[1]-X[2])^3]):

mtaylor~(G(X),[[X[1]=a,X[2]=b]$2],3);

Notice that you need a list of lists as the second argument for the elementwise operation to understand the situation correctly. The problem is that the second argument to mtaylor is a list with as many elements as F(X).

Something like this adapted to your situation:

for k to 5 do FileTools[Rename]( cat("F:/testfile",k,".bla" ), cat("F:/Newtestfile",k,".bla" ) ) end do;

There is a procedure Rename in the FileTools package.

What happens in the call F(2) is that x = 2 is passed to int, so that it reads -int(p(2),2), which is not acceptable since the second argument to int cannot be a number.

You can do like this:

p := x->10*x-20*piecewise(x < 2, 0, x-2);
F1:=unapply(-int(p(x),x),x);

#or like this:

F2:=x->-int(p(t),t=0..x);

# a third version

simplify(-int(p(t),t=0..x));
F3:=unapply(%,x);

The latter two versions use a definite integral, which also ensures that F(0) = 0 (if that is what you want).

For a start you could do something like

with(Student:-NumericalAnalysis);
Lx:=[$-10..10];
Ly:=[i^2$i=-10..10];
for i to nops(Lx) do f(Lx[i]):=Ly[i] end do:
op(4,eval(f));
Quadrature(f(x),x=0..9,partition=9,method=trapezoid);

You will need to do some interpolation to get the integral from 0.5 to 9.1.

With Maple 6 came the LinearAlgebra package. It replaced the old linalg package. Not quite true of course since the linalg package is still part of Maple because of backwards compatibility.

At the same time 'Matrix' and 'Array' were introduced and replacing 'matrix' and 'array', again not quite true.

For 'matrices' the multiplication operator is `&*`, where it is `.`for 'Matrices`.

Try this

A:=LinearAlgebra:-RandomMatrix(2);
B:=LinearAlgebra:-RandomMatrix(2);
A&*B;
C:=evalm(%);
whattype(C);
type(C,Matrix);
type(C,matrix);

For 'Matrices' I have been using `&.`as a neutral operator which is sometimes convenient in exhibiting structures. Then I have extended 'value' to do what eval does below:

A&.B;
eval(%,`&.`=`.`);

restart;
vars:=[a,b,c,d,e,f,g,h,i,j,k]:
for ii from 1 to 10 do
p[ii]:=randpoly(vars,degree=1,homogeneous);
end do:
res:=table():
for var in vars do res[var]:=nops(select(has,[seq(p[ii],ii=1..10)],var)) end do:
eval(res);

I tried this and it worked:

restart;
fd:=["F:/file1.m", "F:/file2.m"];

for i from 1 to nops(fd)
do
cat(var,i):=12345*i;
save cat(var,i),fd[i];
end do;

restart;
fd:=["F:/file1.m", "F:/file2.m"];
for i from 1 to nops(fd)
do
read fd[i];
end do;
seq(cat(var,i),i=1..nops(fd)): %;

Change BCs := [x(0), x(0.2e-1)-0.15e-2] to

BCs := x(0), x(0.2e-1)-0.15e-2;

But you still have problems, but they are not syntactical anymore.

method=bvp[midrich] could be tried, but now you have other problems.

Try

sol := dsolve({BCs, eq}, {x(t)}, type = numeric,method=bvp[midrich],initmesh=1024);

plots:-odeplot(sol,[t,x(t)],0..0.02);

In the standard interface you could use Vector ( = Vector[column]):

Vector([x^2-1,x^3+1]):
factor~(%);