Try this:

restart;
n := 3;
y[0] := x-> 1+x+(1/2)*x^2;
for m from 0 to n do
  y[m+1]:=unapply( y[m](x)+int(-(1/2)*(s-x)^2*(diff(y[m](s),s,s,s)-y[m](s)), s = 0 .. x),x)
end do;
ode:=diff(Y(x),x,x,x)-Y(x)=0;
dsolve({ode,Y(0)=1,D(Y)(0)=1,(D@@2)(Y)(0)=1});
sort(collect(y[n](x),x),x,ascending);
mtaylor(exp(x),x=0,11);
%%-%;

In early versions of Maple all packages were tables. That is true of plots still in Maple 8.

So replace plots:-pointplot with plots[pointplot] . That can even be done in later versions.

Before getting to the parameter you are talking about (actually your worksheet doesn't contain any) you should address the problem that you have: an error message from `dsolve/numeric/bvp/convertsys`.

A solution to this problem is to differentiate the equation not containing the 4th derivatives. In addition use the undifferentiated version to extract an extra boundary condition, because you only have 10, you need 4+4+3 = 11.

restart;
##I should excuse the unnecessarily complicated look of dsys3. I converted to 1-D math input and didn't bother to clean up the mess:
dsys3 := {8*(diff(f2(x), x, x, x, x))+9*(diff(f2(x), x, x))+10*f2(x)+11*(diff(f1(x), x, x, x))+12*(diff(f1(x), x))+13*(diff(f3(x), x, x))+14*f3(x)+f3(x)*f3(x)+(diff(f3(x), x))*(diff(f3(x), x))+(diff(f3(x), x, x))*f3(x) = 0, 16*(diff(f3(x), x, x, x, x))+18*(diff(f3(x), x, x))+19*(diff(f3(x), x, x))+22*(diff(f1(x), x))+23*(diff(f1(x), x))+24*(diff(f2(x), x, x))+25*f2(x)+26*f2(x)+27*f3(x)+29*f3(x) = 0, 2*(diff(f1(x), x, x))+3*(diff(f2(x), x, x, x))+4*(diff(f2(x), x))+6*(diff(f3(x), x))+7*f1(x)+(diff(f3(x), x, x))*(diff(f3(x), x))+(diff(f3(x), x))*f3(x)+(diff(f3(x), x))*f3(x) = 0, f1(0) = 0, f1(1) = 0, f2(0) = 0, f2(1) = 0, f3(0) = 0, f3(1) = 0, ((D@@1)(f2))(0) = 0, ((D@@1)(f2))(1) = 0, ((D@@1)(f3))(0) = 0, ((D@@1)(f3))(1) = 0};
dsys3;
sys,bcs:=selectremove(has,dsys3,diff);
solve(sys,{diff(f3(x),x$4),diff(f2(x),x$4),diff(f1(x),x$3)});
sys[1];
sys[2];#In my case the one I want to differentiate
sys[3];
sys2:=diff(sys[2],x);
#Check that we can solve for the highest derivatives (and defining the new system:
SYS:=solve({sys[1],sys2,sys[3]},{diff(f3(x),x$4),diff(f2(x),x$4),diff(f1(x),x$3)});
bcs;
#Finding one more boundary condition:
convert(sys[2],D);
eval(%,x=0);
bcs1:=eval(%,bcs);
## Now we can solve the problem:
dsol:= dsolve(SYS union bcs union {bcs1}, numeric);
plots:-odeplot(dsol,[x,f3(x)],0..1);
#All 3 functions are dead zero: Not surprising with the homogeneous boundary conditions.

So the parameter you are talking enters where?
May we assume that you want a nontrivial solution?
You could experiment with giving an approximate solution as in this rather crude attempt:
dsol:= dsolve(SYS union bcs union {bcs1}, numeric,
       approxsoln=[f1(x)=100*x*(1-x),f2(x)=100*sin(Pi*x),f3(x)=100*x*(1-x)]);
plots:-odeplot(dsol,[[x,f1(x)],[x,f2(x)],[x,f3(x)]],0..1);

The particular output here I wouldn't trust right away.

With the assumptions x>0 and t>0 we can do:
restart;
sys := [v*diff(u(x,t), x) + diff(u(x,t), t) = 0, u(x,0) = exp(-x), u(0,t) = sin(t)];
res:=pdsolve(sys[1]);
eval(res,{t=0,u(x,t)=exp(-x)});
eval(res,{x=0,u(x,t)=sin(t)});
#Thus we get:
_F1:=y->piecewise(y<0,exp(v*y),sin(y));
#So that the solution is
res;
sol1:=collect(convert(res,Heaviside),Heaviside);
#sol1 is the same as
sol:=u(x, t) = exp(t*v-x)+Heaviside(t-x/v)*(sin(t-x/v)-exp(t*v-x));
simplify(sol-sol1) assuming t>x/v;
simplify(sol-sol1) assuming t<x/v;

I think I'm beginning to understand your loop. I shall use the following version. I use 'Linv' I defined in my other answer. I removed the option 'continuous' as that was basically put there to make Linv work with symbolic y. But you use it on concrete values for y.
I take it that in the line defining R[m] with g(y) you actually mean g(y[m-1]). That is what I use below anyway.
I only go to n=3 since only y[1] and y[2] are actually computed all the way, y[3] is still a double integral, which, however can be evaluated numerically.
I compare the results to solving numerically by using dsolve/numeric.
The result for y[3] is not bad.

restart;
Linv:=z->unapply(int(s^(-2)*int(t^2*z(t),t=0..s),s=0..x),x);
Y[0]:=0:
g:=y->exp(-y):
n:=3: ##
for m from 1 to n do
  YF:=unapply(Y[m-1],x);
  R:=Y[m-1]+Linv(g@YF)(x); ##This line and the one below could be shortened
  Y[m]:=Y[m-1]-R               ## Actually Y[m] is just defined as -Linv(g@YF)(x);
end do;
eval(Y[3],x=0.5); #Test computation
plot(Y[3],x=0..1,adaptive=false,numpoints=25);#Takes a little while
plot([Y[1],Y[2]],x=0..1);#Fast
# Numerical solution as a comparison: I consider this basically as the "exact" solution.
ode:=diff(y(x),x,x)+2/x*diff(y(x),x)+exp(-y(x))=0;
res:=dsolve({ode,y(0)=0,D(y)(0)=0},numeric,output=listprocedure);
Ynum:=subs(res,y(x));
plots:-odeplot(res,[x,y(x)],0..1); #The "exact" solution
plot(Ynum(x)-Y[3],x=0..1,adaptive=false,numpoints=25,caption="Error for Y[3]"); Y3err:=%:
plots:-display(Array([seq(plots:-odeplot(res,[x,y(x)-Y[i]],0..1,caption=typeset("Error for ",'Y[i]')),i=1..2),Y3err]));
##########Simpler version ##############
##Using that after changing the order of integration in int( 1/s^2*int(t^2*z(t),t=0..s), s=0..x) we can perform the s-integration to obtain int(t^2*(1/t-1/x)*z(t),t=0..x).
#So we get for n = 3:
Y[0]:=0:
for m from 1 to 3 do
  YF:=unapply(Y[m-1],x);
  Y[m]:=-int(t^2*(1/t-1/x)*exp(-YF(t)),t=0..x)
end do;

There may be a smarter way, but this works on my example:

f:=1.23456789*10^7/(9.87654321*10^6*a+6.999988888*10^8*b);
evalindets(f,float,x->x*10^(-6));
##I'm dividing all floats by 10^6.

You can implement the operator L this way:

restart;
L:=y->unapply(convert(x^(-2)*diff(x^2*diff(y(x),x),x),D),x);
# Thus L(y) is a function:
L(y)(t);
L(y)(4);
## To find an implementation of an inverse, we first try stepwise:
t^2*L(y)(t);
int(convert(%,diff),t=0..s,continuous);
int(s^(-2)*%,s=0..x,continuous);
## Notice that we must assume that the function space on which L is defined satisfies y(0)=0 for all y.
##In one line
int(s^(-2)*int(t^2*L(y)(t),t=0..s,continuous),s=0..x,continuous);
convert(%,diff);
## Thus Linv is (under the assumption made):
Linv:=z->unapply(int(s^(-2)*int(t^2*z(t),t=0..s,continuous),s=0..x,continuous),x);
(Linv@L)(y)(x);
convert(%,diff);
(L@Linv)(y)(x);
##############
As for the loop I don't know what you are trying to do.

You can use the 'known' option:

restart;
ff:= x -> if not type(x,numeric) then 'procname(_passed)' else evalf(Int(exp(t), t= 0 .. x, method = _d01ajc) + 1) end if;
ff(x); #test
ff(1); #test
ode:=D(f)(x) = ff(x);
sol:=dsolve({ode, f(0)=1}, numeric, known=ff,method=gear);
sol(1); #Test

You can replace the method 'gear' by lsode, dverk78, or by several classical, e.g. classical[rk4].

####### Notice the need for returning unevaluated.
Although for a name x
evalf(Int(exp(t), t= 0 .. x, method = _d01ajc) + 1);
returns an unevaluated Int-integral, this integral has been stripped of its option (and evalf is gone).

As Rouben points out the first equation is a pde, so use pdsolve

pde := a*(diff(M(a, b, c), a))+a*b*(diff(M(a, b, c), b))*c = 0;  # ..=0 is assumed, but better put it there
pdsolve(pde);

In the result  M(a, b, c) = _F1(c, b*exp(-c*a)),  _F1 is an arbitrary function of two variables, thus e.g.

_F1:= (x,y) -> x^2*sin(y);
would be OK. Then the solution would look like this
M(a, b, c) = c^2*sin(b*exp(-c*a))

#######  The second equation has two functions differentiated, a(t) and b(t). But you only got that one equation.
Thus you must decide which is the unknown function.
If it is a(t) then do:

ode := a(t)*(diff(a(t), t))+a(t)*b(t)*(diff(b(t), t))*c(t) = 0;
dsolve(ode,a(t));

Spurred by acer's workaround and trying to strip `convert/global` to the essential part I found a very simple workaround. 
It appears that no locals occur at all in the examples we have considered, so that part could be stripped. However, there were several uses of indets inside the procedure.
I ended up with the following procedure to be used instead of `convert/global` in the present context.

cg4:=proc(nom) local a; indets(a,symbol); nom end proc;

It is just the identity map. However, it must have a side effect that is useful. I made 'a'  local, but that is not essential. The type 'symbol' might as well have been 'float' or 'function', or whatever (I haven't checked more).

Using cg4 in Carl's example:
indets(expr, HODD &under (cg4@`convert/D`));
gives the expected result.
#### Here are several of the examples mentioned in my comments above. All are successfully handled by cg4.
It even does better than `convert/global` in the last example w3.

restart;
cg4:=proc(nom) local a; indets(a,symbol); nom end proc;
w1:=diff(u(x),x)+D(u)(x)+Diff(u(x),x);
indets(w1, specfunc(D(u)) &under (cg4@`convert/D`));
indets(w1,specfunc(diff) &under (cg4@`convert/diff`));
indets(w1,specfunc(Diff) &under (cg4@`convert/Diff`));
########################################
w2:=[int(f(x),x),Int(f(x),x)];
indets(w2,specfunc(int) &under (cg4@`convert/int`));
indets(w2,specfunc(int) &under (cg4@value));
indets(w2,specfunc(Int) &under (cg4@`convert/Int`));
#########################################
w3:=[f,%f];
indets(w3,identical(f) &under (cg4@value));
indets(w3(x),specfunc(f) &under (cg4@value));
indets(w3,identical(f) &under (`convert/global`@value)); #Only %f
indets(w3,identical(%f) &under (`convert/global`@curry(subs,f=%f))); # Only f
indets(w3,identical(%f) &under (cg4@curry(subs,f=%f))); #OK!
indets(w3(x),specfunc(%f) &under (`convert/global`@curry(subs,f=%f))); #OK!
indets(w3(x),specfunc(%f) &under (cg4@curry(subs,f=%f))); #OK!
##########################################
## Finally an oddity probably not unrelated to why cg4 works:
indets(w1, specfunc(D(u)) &under (`convert/D`@@2)); ### Works!
indets(w1, specfunc(diff) &under (`convert/diff`@@2)); # Only two of them
indets(w1, specfunc(Diff) &under (`convert/Diff`@@2)); # Only two of them

Use combine to obtain what you want:
combine(1-2*sin(x)^2);

simplify has a general preference for cos over sin. That doesn't mean however, that it turns sin into cos at all costs:

simplify(sin(x));
##Try also
simplify(1-2*sin(x)^2,size);

simplify doesn't necessarily get you the simplest result in the common sense of the word 'simplify'. Try as another example

expand((x+y)^3);
simplify(%);
factor(%);

You can use eval[recurse]:

restart;
f:=x -> 1+x^2/n;
g:=f@@n;
g(x);
eval[recurse](g(x),{n=3});

In eval(g(x),n=3) what happens is that g(x), which just evaluates to (f@@n)(x), is evaluated to (f@@3)(x) and the result will still contain the symbol n.

You have two different e's in the very first image. In the (1,1) element one 'e' is due to the use of e^(-t), which should have been exp(-t), the other 'e' is printed slightly different and is due to the correct use in exp(-4*t).  Similarly for the other 3 matrix elements.

Here is another approach at solving the integral equation.
In its present version it is rather crude. I have to admit to knowing nothing about numerical solutions to integral equations. But here it is.
It assumes that what you want to solve is actually the equation referred to as eq just below. Furthermore that you want y(0) = 1.
##The first version was flawed and has been replaced by the following:
restart;
Digits:=15:
macro(AInt=Student:-Calculus1:-ApproximateInt);
eq:=y(t)=1-h*Int(JacobiTheta3(0,exp(-Pi^2*s))*y(t-s)^4,s=0..t);
eq2:=subs(Int=rcurry(AInt,partition=25),eq); #Default for AInt is method=midpoint
#eq2:=y(t)=1-h*AInt(JacobiTheta3(0,exp(-Pi^2*s))*y(t-s)^4,s=0..t,method=midpoint,partition=25):#Same as above
h := 0.65e-4;
n:=10:
Y:=t->1+add(c[i]*t^i,i=1..n);
eqt:=eval(eq2,y=Y):
T:=[seq(0.1..1,0.1)];
nops(T); #Same as n
eqs:={seq(evalf(eval(IntegrationTools:-Expand(eqt),t=t0)),t0=T)}:
sol:=fsolve(eqs);
eval((rhs-lhs)~(eqs),sol); #Checking collocation
Ymono:=eval(Y(t),sol); #The solution using the monomials t^i
plot(Ymono,t=0..1,size=[1800,400]);

I tried your approach with infinity (inf, in your case 100) replaced by 3. The loop went from just 1 to N with N=2.
You can imagine the explosion in the number of terms you have to deal with if inf and/or N are chosen much larger.

restart;
h := 0.65e-4;
inf:=3:
Theta3:=t->add(exp(-m^2*Pi^2*t),m=-inf..inf);
y[0]:=t->1;
N:=2:
for i from 1 to N do
  y[i]:=unapply(combine(expand(1-h*int(Theta3(s)*y[i-1](t-s)^4,s=0..t))),t)
end do;
plot(y[1](t)-y[2](t),t=0..10);
plot(y[2](t),t=0..10);

Briefly I also looked at a change of variable with x = exp(-t*Pi^2) thereby transforming the t-interval 0..infinity into the x-interval 0..1. That doesn't get by the explosion of terms of course.