That's strange.  I don't get that result after Delta := 3 either:

> filtre:=proc(g) local t;
   return x->(1/(2*Delta)*int(g(t),t=(x-Delta)..(x+Delta)));
   end;
  Delta := 3;
  filtre(Dirac)(x);

1/6*Heaviside(x+3)-1/6*Heaviside(x-3)

What version of Maple are you using?  Can you upload a worksheet in which you get your result?

That's strange.  I don't get that result after Delta := 3 either:

> filtre:=proc(g) local t;
   return x->(1/(2*Delta)*int(g(t),t=(x-Delta)..(x+Delta)));
   end;
  Delta := 3;
  filtre(Dirac)(x);

1/6*Heaviside(x+3)-1/6*Heaviside(x-3)

What version of Maple are you using?  Can you upload a worksheet in which you get your result?

It works for me in Maple 11: I get:

> filtre:=proc(g) local t;
   return x->(1/(2*Delta)*int(g(t),t=(x-Delta)..(x+Delta)));
  end;
  filtre(Dirac)(x);

1/2*1/Delta*(Heaviside(x+Delta)-Heaviside(x-Delta))

It works for me in Maple 11: I get:

> filtre:=proc(g) local t;
   return x->(1/(2*Delta)*int(g(t),t=(x-Delta)..(x+Delta)));
  end;
  filtre(Dirac)(x);

1/2*1/Delta*(Heaviside(x+Delta)-Heaviside(x-Delta))

Numerical integration is probably not the best application for Pade, but here's what I get.  Note that the fair comparison with pade(..., [4,4]) would by taylor(..., 9), because the Pade approximation is obtained from that Taylor series.

> with(numapprox):
   f:= sec(x)/(x+2);
   T := convert(taylor(f, x=0, 9), polynom); 
   P := pade(f, x, [4,4]);
   evalf([Int(f,x=-1..1), Int(T,x=-1..1), Int(P,x=-1..1)]);

 [1.369112680, 1.364325672, 1.369098559]

So here the Pade is much better than the Taylor.

But for a nicer example to show off the power of Pade, try a function that has poles in the complex plane, and an interval on which the Taylor series won't converge (the radius of convergence being the distance to the nearest pole).

> f := x/(exp(x) + exp(-x));
   T:= convert(taylor(f, x=0, 19), polynom);
   P:= pade(f, x, [9,9]);
   evalf([Int(f, x=0..5), Int(T, x=0..5), Int(P,x=0..5)]);

[.8755384560, 88407513.49, .8755387258]

Numerical integration is probably not the best application for Pade, but here's what I get.  Note that the fair comparison with pade(..., [4,4]) would by taylor(..., 9), because the Pade approximation is obtained from that Taylor series.

> with(numapprox):
   f:= sec(x)/(x+2);
   T := convert(taylor(f, x=0, 9), polynom); 
   P := pade(f, x, [4,4]);
   evalf([Int(f,x=-1..1), Int(T,x=-1..1), Int(P,x=-1..1)]);

 [1.369112680, 1.364325672, 1.369098559]

So here the Pade is much better than the Taylor.

But for a nicer example to show off the power of Pade, try a function that has poles in the complex plane, and an interval on which the Taylor series won't converge (the radius of convergence being the distance to the nearest pole).

> f := x/(exp(x) + exp(-x));
   T:= convert(taylor(f, x=0, 19), polynom);
   P:= pade(f, x, [9,9]);
   evalf([Int(f, x=0..5), Int(T, x=0..5), Int(P,x=0..5)]);

[.8755384560, 88407513.49, .8755387258]

By "sigma format", I assume you mean something like this:

> Slode[mhypergeom_formal_sol](((D@@2)(y))(x)-x*y(x) = 0, y(x));

_C[0]*GAMMA(2/3)*Sum(9^(-n)/GAMMA(n+1)/GAMMA(n+2/3)*x^(3*n),n = 0 .. infinity)+2*_C[1]*Pi/GAMMA(2/3)*Sum(3^(-3/2-2*n)/GAMMA(n+4/3)/GAMMA(n+1)*x^(3*n+1),n = 0 .. infinity)

As for more terms in the result of dsolve(...,series), you can set the environment variable Order.  For example:

> Order := 30:
  dsolve(((D@@2)(y))(x)-x*y(x) = 0,y(x),series);

y(x) = y(0)+D(y)(0)*x+1/6*y(0)*x^3+1/12*D(y)(0)*x^4+1/180*y(0)*x^6+1/504*D(y)(0)*x^7+1/12960*y(0)*x^9+1/45360*D(y)(0)*x^10+1/1710720*y(0)*x^12+1/7076160*D(y)(0)*x^13+1/359251200*y(0)*x^15+1/1698278400*D(y)(0)*x^16+1/109930867200*y(0)*x^18+1/580811212800*D(y)(0)*x^19+1/46170964224000*y(0)*x^21+1/268334780313600*D(y)(0)*x^22+1/25486372251648000*y(0)*x^24+1/161000868188160000*D(y)(0)*x^25+1/17891433320656896000*y(0)*x^27+1/121716656350248960000*D(y)(0)*x^28+O(x^30)

By "sigma format", I assume you mean something like this:

> Slode[mhypergeom_formal_sol](((D@@2)(y))(x)-x*y(x) = 0, y(x));

_C[0]*GAMMA(2/3)*Sum(9^(-n)/GAMMA(n+1)/GAMMA(n+2/3)*x^(3*n),n = 0 .. infinity)+2*_C[1]*Pi/GAMMA(2/3)*Sum(3^(-3/2-2*n)/GAMMA(n+4/3)/GAMMA(n+1)*x^(3*n+1),n = 0 .. infinity)

As for more terms in the result of dsolve(...,series), you can set the environment variable Order.  For example:

> Order := 30:
  dsolve(((D@@2)(y))(x)-x*y(x) = 0,y(x),series);

y(x) = y(0)+D(y)(0)*x+1/6*y(0)*x^3+1/12*D(y)(0)*x^4+1/180*y(0)*x^6+1/504*D(y)(0)*x^7+1/12960*y(0)*x^9+1/45360*D(y)(0)*x^10+1/1710720*y(0)*x^12+1/7076160*D(y)(0)*x^13+1/359251200*y(0)*x^15+1/1698278400*D(y)(0)*x^16+1/109930867200*y(0)*x^18+1/580811212800*D(y)(0)*x^19+1/46170964224000*y(0)*x^21+1/268334780313600*D(y)(0)*x^22+1/25486372251648000*y(0)*x^24+1/161000868188160000*D(y)(0)*x^25+1/17891433320656896000*y(0)*x^27+1/121716656350248960000*D(y)(0)*x^28+O(x^30)

This is not a problem I've encountered with any version of Maple. 
What operating system is this under?
Are you sure the directory you're saving it to is one where you have permission to write files?

That has a logarithmic singularity at x = cos(2), so it can't be an antiderivative of
arccos(x)/(1+x^4).

To avoid those commands that Maple seems to find difficult, we can solve the nullcline equations explicitly for v.

This was dv/dt:

> vp := -1/200000*w-7/100*v-3/1000+2000/w*v^2+80/w*v;

> Vv := solve(vp, v);

Vv := 7/400000*w-1/50+1/400000*(449*w^2+128000*w+64000000)^(1/2), 7/400000*w-1/50-1/400000*(449*w^2+128000*w+64000000)^(1/2)

The first of these two is the one we are most interested in here, which is 0 at w=0.

This was dw/dt:

> wp := 1/10*w-v*w-2000*v;

> Vw := solve(wp, v);

Vw := 1/10*w/(w+2000)

It is easy to check that Vv[1] and Vw are increasing for w > 0, with
Vv[1] < Vw for 0 < w < 171.688309 approximately.  For w > 0 and
Vw[2] < v < Vv[1], v' < 0 and w' > 0.  For v< Wv[1] < v < Vw, v' > 0 and w' > 0.
For v > Vw, v' > 0 and w' < 0.  So the picture really is as I described it.

You're looking at this system:

> sys:=convert([xdot,ydot,zdot],rational);

sys := [diff(x(t),t) = 20*x(t)^(1/2)-y(t), diff(y(t),t) = (10/x(t)^(1/2)-z(t))*y(t), diff(z(t),t) = (z(t)-1/25)*(1/20-z(t))+(z(t)-3/50)*(10/x(t)^(1/2)-z(t))*y(t)/(20*x(t)^(1/2)-y(t))]

Because of the 20*sqrt(x(t))-y(t) in the denominator, there won't be an actual equilibrium point, and indeed you have singularities where the differential equation becomes undefined.  But you're interested in seeing what happens near such points. 

I think it's helpful to use a change of variables: u = sqrt(x) (to remove those annoying square roots), v = z - 10/sqrt(x) (so dy/dt = 0 when v = 0), w = y - 20 sqrt(x) (so the singularity is along w = 0).
 

> PDEtools[dchange]({x(t)=u(t)^2,y(t)=w(t)+20*u(t),z(t)=v(t)+10/u(t)},
      sys,[u(t),v(t),w(t)]);
 simplify(%) assuming u(t)>0;  
  newsys:= expand(solve(%, 
       {diff(u(t),t), diff(v(t), t), diff(w(t),t)})) ;

newsys := {diff(u(t),t) = -1/2*w(t)/u(t), diff(w(t),t) = 10*w(t)/u(t)-v(t)*w(t)-20*v(t)*u(t), diff(v(t),t) = -5/u(t)^3*w(t)+3/100*v(t)-10/u(t)*v(t)+9/10*1/u(t)-100/u(t)^2-1/500+20*u(t)/w(t)*v(t)^2+200/w(t)*v(t)-6/5*u(t)/w(t)*v(t)}

When w and v are small (but u is not), u will change relatively slowly.  So maybe it will be helpful to see the phase plane of v and w with u kept at a constant value, say 100.

> eval(newsys, u(t) = 100);
   sys100:= select(has,%, diff);
   with(DEtools): DEplot(sys100,[w(t),v(t)],t=-1..1,
      w=-0.1 .. 0.1,v=-0.1..0.1,arrows=medium);

This looks rather strange, because the arrows seem to do a 180 degree turn as they cross the w axis.  But here's a closer look, with the isoclines added.

> vp := subs(sys100, v(t)=v, w(t)=w, diff(v(t),t));
   wp := subs(sys100, v(t)=v, w(t)=w, diff(w(t),t));
  Digits:= 15:
  with(plots):
  display([
     implicitplot([vp,wp], w=-0.1  .. 0.1, v = -0.001 .. 0.001,
         colour=[red,blue],gridrefine=3),
     DEplot(sys100,[w(t),v(t)],t=-1..1,w=-0.1 .. 0.1,v=-6e-6..6e-6,
          arrows=medium)]);

The wedge-shaped region between the red and blue isoclines on the right
is an "antifunnel" in John Hubbard's terminology: trajectories can leave it but
can't enter it, because the blue isocline is a lower fence and the red isocline is
an upper fence.  Trajectories come out of the origin to the right inside the
antifunnel, and can leave it through the red isocline (then moving off to the right with v negative) or through the blue isocline (then moving up, w decreasing but staying positive).  On the left, trajectories coming up from the left, below the blue isocline,  leave the picture going down and to the left, just above the red isocline.

This is not the usual setup for a phase portrait, because it is only a first-order
differential equation: kdot = diff(k(t),t) is a function of k, so you'd just have a curve
in the k - kdot plane.  For a phase plane you'd want a system of two first-order
DE's, or a single second-order DE which would be converted into such a system.
You could do something like this, I guess:
 

> with(plots):
  f := k -> 0.7 * k^0.3 - 0.13*k;
  P1:= plot(f(k), k = 0 .. 15):
  P2 := seq(arrow([j, f(j)], [f(j),0]), j=0..15):
  display([P1,P2],labels=[k,k(dot)]);

This is not the usual setup for a phase portrait, because it is only a first-order
differential equation: kdot = diff(k(t),t) is a function of k, so you'd just have a curve
in the k - kdot plane.  For a phase plane you'd want a system of two first-order
DE's, or a single second-order DE which would be converted into such a system.
You could do something like this, I guess:
 

> with(plots):
  f := k -> 0.7 * k^0.3 - 0.13*k;
  P1:= plot(f(k), k = 0 .. 15):
  P2 := seq(arrow([j, f(j)], [f(j),0]), j=0..15):
  display([P1,P2],labels=[k,k(dot)]);

Interesting: the expected number of repetitions of length 7 in 10000 random digits is approximately 5, so this is a bit low, though not unusually low.

The problem with a larger string is that Repeats will take a long time and return a huge sequence, with about N^2/20 entries for a random string of N digits.  That's already 5 million for your 10000 digits.

Perhaps a better idea would be this.  Suppose you're interested in repetitions of length m in a string S of N digits.   Let T be a table indexed by m-digit strings.  For each m-digit string w in S, we check if T[w] exists; if it does, put a new entry in the table "mreps", otherwise make a new entry T[w] = (index where w starts in S).
 

So here is how I'd find repetitions of length 9 in the first 10^5 digits of Pi.  Again the expected number of these is about 5.

> with(StringTools):  
   PiDigits:= Delete(convert(evalf(Pi,10^5),string),2..2):
  for ind from 1 to 10^5-9 do
   w:= PiDigits[ind..ind+8];
   if assigned(T[w]) then mreps[w]:=[T[w],ind]
     else T[w]:= ind
   fi
 od: 
 eval(mreps);

TABLE(["041021944" = [21762, 75871], "134926158" = [69598, 86283], "201890888" = [30628, 33845]])