I assume you're using the Student[NumericalAnalysis] package.  The option in Euler that controls the step size is numsteps, which gives the number of steps used.  If you're going from t = a to t = b in n steps, the step size will be (b-a)/n.  Thus for step size h = 0.05 from a = 0 to b = 0.4, you would use

Euler(DE2, y(0) = 1, t = .4, numsteps = 8);

Is your curve defined by an expression or function?  You should be able to take an interval around each maximum, small enough to contain only one critical point, and use fsolve to find the critical point in that interval.  Something like:

seq(fsolve(diff(F, x) = 0, x = J), J = [ 687 .. 690, 705 .. 710, 736 .. 740, 758 .. 760 ]);

> q:= (-1+du1(t))/(C*L);
   sort(-1+du1(t), du1(t));
  
And then

> q;

You might try quickplot3d from my Maple Advisor Database, <http://www.math.ubc.ca/~israel/advisor/>.

For example, after installing the database you should be able to do this:

M := Matrix(readdata("sample.dat",float,3), datatype=float[8]);
quickplot3d(M, style=point, symbol=cross);

You may find the Classic GUI better than Standard at handling plots with large numbers of points.

There are two distinct types in Maple: matrix and Matrix.  matrix is used in the old linalg package, and uses evalm(A &* B) for multiplication, Matrix is used in LinearAlgebra, and uses A . B.  The old versions are still included for the sake of backward compatibility, but I highly recommend using the new versions, Matrix and A . B. 

1) I presume you mean angledata rather than thetadata
2) missing ] in definition of plotseq
3) plots:-display rather than just display
4) your pendulum seems to be upside-down
5) it may work out ok in this case, but in general there's no guarantee that sol(i/30)[2] will be the equation you want.  A better style would be to use subs(sol(i/30), theta(t)) instead.

Why not use an indexed variable instead?

assume(a,real);
assume(b,real):

parameters:=[a,b];

i:=0;
for what in parameters do
temp:=i;
crib[what]:= temp;
i:=i+1;
end do;

crib[a], crib[b];

    0, 1

First you have to write your expression correctly.  Don't use square brackets for grouping: those are for lists and indices.  And use * for multiplication.  Also, it's usually a bad idea to use the same name in both indexed and un-indexed forms (x and x[n]), so I'll replace x[n] by X[n]. 

> y(t):= (1/3*h^2)*((-x+X[n]+h)*(-x+X[n]+3*h))*y[n] + (1/3*h^2)*((-X[n]+x)*(-x+X[n]+4*h))*y[n+1]-(1/3*h)*((-X[n]+x)*(-x+X[n]+3*h))*f[n+2];

> convert(series(y(t), x = X[n]),polynom);

Actually implicitplot may not be the best way to produce a solid-colour region.  Why not just use plottools:-rectangle?

a:=plot(x^3-4*x,x=-3..3):
b:= plottools:-rectangle([-5,-5],[5,5],colour=black):
plots:-display(a,b,view=[-5..5,-5..5],axis=[colour=green]);

The point, I think, is that the equation y = 1 + A t + B t^2 + C t^3 can be rewritten as

t = (y - 1)/(A + B t + C t^2)

so you can try finding approximate solutions by iterating

t[n] = (y - 1)/(A + B t[n-1] + C t[n-1]^2)

However, I don't think it's reasonable to expect a "closed-form" solution to this recurrence, and that's what rsolve would look for.  Even in rather simple cases, it fails to find any.

> rsolve({t(n)= 1/(1  + t(n-1)^2),t(1) = 1}, t(n));

(returns unevaluated)

 So forget rsolve, take an example, and just iterate a few times to see if it appears to converge (I suggest using floating point)

The problem is the semicolon at the end of the first line.  Remove it and it should parse OK.
However, your code has other problems.  For one thing, you're evaluating everything at
(tinit, xinit, yinit) instead of (t[i-1],x[i-1],y[i-1]).

As hirnyk suggests, it may be useful to take the logarithms of your y data.  However, I might suggest a quartic polynomial fit for log(y) as a function of x. 

  > datax := [1, 1.019, 1.040, 1.059, 1.079, 1.100, 1.120, 1.140, 1.160, 1.180, 1.200];
     datay := [2.626, 3.007, 2.955, 2.982, 2.940, 3.527, 6.433, 16.172, 36.477, 106.440, 394.928];
     logdatay:= map(ln, datay);
     F:=Statistics[Fit](add(a[i]*x^i,i=0..4),datax,logdatay,x);

F := -15286.3322412965+55513.2171569571*x-75334.3325224527*x^2+45263.1842157640*x^3-10154.7703713389*x^4

Thus the suggestion is that y = exp(F) with F as above.

 > with(plots): display(plot(F,x=1..1.2),pointplot(zip(`[]`,datax,logdatay)));      

And here's a version of the procedure that plots the connected component containing 1+i in black.

nb:= proc(r) local i,j,nb1; nb1:={seq(seq(i+I*j,j=0..floor(sqrt(r^2-i^2))),i=0..r)};
  map(t -> (t,I*t,-t,-I*t),nb1 minus {0}) end proc;;
    myF:= proc(R,r)
       uses GraphTheory;
       local nbd,G,E,Gs,Es,CC1, G1, E1, E2, Gr, CC;
       G:= map(g -> (g,I*g,-g,-I*g),
               select(GaussInt[GIprime],{seq(seq(x+I*y,y=0 .. floor(sqrt(R^2 - x^2))),x=0..R)}));
       nbd:= nb(r);
       Gs:= convert(map(convert,G,string),list);
       E:= `union`(seq(map(t -> {g,t},select(member,map(t -> t+ g,nbd),G)),g=G));
       Es:= map(e -> map(convert,e,string), E);
       Gr:= Graph(Gs,Es);
       CC:= ConnectedComponents(Gr);
       G1:= map(parse,op(select(has,CC,"1+I")));
       E1,E2:= selectremove(e -> member(e[1],G1),E);
       E1:= map(e -> [[Re(e[1]),Im(e[1])],[Re(e[2]),Im(e[2])]], E1);
       E2:= map(e -> [[Re(e[1]),Im(e[1])],[Re(e[2]),Im(e[2])]], E2);
       plots[display](PLOT(CURVES(op(E1),COLOUR(RGB,0,0,0))),
                      PLOT(CURVES(op(E2),COLOUR(RGB,1,0,0))),
 plottools[circle]([0,0],R,colour=blue),axes=none);
    end proc;

> myF(100,3);

Here's my try:

> nb:= proc(r) local i,j,nb1; nb1:={seq(seq(i+I*j,j=0..floor(sqrt(r^2-i^2))),i=0..r)};
 map(t -> (t,I*t,-t,-I*t),nb1 minus {0}) end proc;;
   myF:= proc(R,r) 
      local G, nbd, E;
      G:= map(g -> (g,I*g,-g,-I*g),
              select(GaussInt[GIprime],{seq(seq(x+I*y,y=0 .. floor(sqrt(R^2 - x^2))),x=0..R)}));
      nbd:= nb(r);
      E:= `union`(seq(map(t -> [[Re(g),Im(g)],[Re(t),Im(t)]],select(member,map(t -> t+ g,nbd),G)),g=G));
      plots[display](PLOT(CURVES(op(E),COLOUR(RGB,1,0,0))),plottools[circle]([0,0],R,colour=blue),axes=none);
   end proc;

 myF(100, 3);

If we write diff(y(x),x) = v, then v^2 - a*y^2 + y^4/2 - y^6/3 is constant on solutions.  Here are some trajectories in the y-v plane.  Note the separatrices in blue: solution curves going into or out of the critical points (y,v)  = (b,0) and (-b,0), where b = 1.311737612 is a real root of y^4 - y^2 - a.

> with(DEtools): a:= -1.24:
   sys:= {diff(y(x),x) = v(x), diff(v(x),x) = a*y(x) - y(x)^3 + y(x)^5};
   DEplot(sys,[y(x),v(x)],x=-5..5,
      [[y(0)=1.31,v(0)=0],[y(0)=-1.31,v(0)=0],[y(0)=1.32,v(0)=0],[y(0)=-1.32,v(0)=0], 
       [y(0)=0.9,v(0)=0],[y(0)=2,v(0)=0],[y(0)=-2,v(0)=0]],y=-3..3,v=-4..4,linecolor=[blue$4,green$3]):

Now bcs2 should correspond to solutions that start on the line v = 0 and then hit y = 0 after "time" L.
Note that this can only happen in the region of closed trajectories, where by symmetry L is (n/2+1/4) times the period of the solution for some integer n.  The period of the solution starting at (y,0) (for y from 0 to approximately 1.31) s obtained by the following procedure:

> Y := subs(dsolve({diff(y(x),x$2) - a*y(x) + y(x)^3 - y(x)^5  = 0,y(0)=y0,D(y)(0)=0},y(x),numeric, parameters=[y0],output=listprocedure). y(x));
   P:= proc(yy) Y(parameters=[y0=yy]); forget(fsolve); 4*fsolve(Y(x)=0,x=0..4) end proc;