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~(%);

Here is an attempt that looks like it is working. Maybe there is a better way.

I begin by changing variables. The change x= z*d/2 is just for convenience.

The idea behind the other change y(x)=u(z)*y0 is that y0 should be considered y(0) later. Your problem as stated obviously has the trivial solution y(x) = 0 for all x. You want to find nontrivial solutions. Therefore that change.

The problem then is: Given d find y0 so u(0)=1.

restart;
ode:= diff(y(x),x$2) + 1/x *diff(y(x),x) + 0.6*y(x) + y(x)^3 - y(x)^5  = 0;
PDEtools:-dchange({x=z*d/2,y(x)=u(z)*y0},ode,[z,u]);
ode2:=expand(%*d^2/4/y0);
#New boundary conditions
bcs2:=D(u)(0) = 0,   u(1) = 0.8;
#Don't forget option remember in the following procedure, which outputs u, so if you want y (as a function of z) then you must multiply by y0.
pu:=proc(dd,y00) option remember; local ode3,U,L;
 if not type([dd,y00],list(numeric)) then return 'procname(dd,y00)' end if;
 ode3:=eval(ode2,{d=dd,y0=y00});
 L:=dsolve({ode3,bcs2},numeric,method=bvp[midrich],output=listprocedure);
 subs(L,u(z))
 end proc:

plot(pu(80,1.19),0..1);
seq(pu(80,1.192+c/10000)(0),c=0..9);
plot(1.1925*pu(80,1.1925),0..1);
#I didn't succeed in getting fsolve to produce anything, but see implicitplot below.
fsolve(pu(80,y0)(0)=1,y0=1.9);
plots:-animate(plot,[pu(80,y0),0..1,thickness=2],y0=1.1..1.3,frames=10);
plots:-implicitplot(pu(d,y0)(0)=1,d=2..100,y0=.8..1.2);

#Added after lunch (extracting y0-d-information from implicitplot and using that y0 is constant eventually):

imp:=plots:-implicitplot(pu(d,y0)(0)=1,d=2..10,y0=.8..1.2): imp;
A:=op(indets(imp,Array)):
A[99,1..2];
Y0:=unapply(piecewise(d<=10,CurveFitting:-Spline(A,d),A[99,2]),d):
plot(Y0(d),d=2..100);
plot([seq(Y0(k)*pu(k,Y0(k)),k=[2.1,2.2,2.5,3,5,10,100])],0..1,view=0..1.2);