Example:

f:=1/(sin(x)-1/2);
discont(f,x);

                        /1      2                  \
                       { - Pi + - Pi _B1 + 2 Pi _Z1 }
                        \6      3                  /
Taken from the help page in Maple 12 which happened to be what was within my reach.

10..15; #A range
$10..15; #a sequence
i$i=10..15; #sequence
seq(i,i=10..15); # sequence
seq(10..15); #sequence

For information about $ see
?$

$ is often used in this context:
diff(f(x),x$4);
same as
diff(f(x),x,x,x,x);

Make no assignments, but use `= ` and  set notation.

restart;
N:=3:
system1:={u[0,0]=(1/2)*(u[1,0]+u[0,1]),
u[0,N+1]=(1/2)*(u[0,N]+u[1,N+1]),
u[N+1,N+1]=(1/2)*(u[N+1,N]+u[N,N+1]),
u[N+1,0]=(1/2)*(u[N,0]+u[N+1,1])};
system2:={seq(u[0, j] = (1/4)*(u[1,j]+u[1,j]+u[0, j-1]+u[0,j+1]-f[0,j]*h^2),j=1..N)};
system3:={seq(u[N+1, j] = (1/4)*(u[N, j]+u[N-1,j]+u[N+1, j-1]+u[N+1, j+1]-f[N+1,j]*h^2),j=1..N)};
 eqs := { seq(seq(Stencil[1](h,i,j,u,f),i=1..N),j=1..N)};
sys:=`union`(eqs,system3,system2,system1);

Although you give some context it is not very specific, at least not for someone (like me) who doesn't know anything about quantum tunneling.
Just in case this has to do with solving odes I give this example (which is taken out of my head):

sys:={diff(psi[1](x),x,x)=psi[2](x),diff(psi[2](x),x)=psi[1](x)+psi[2](x),psi[1](0)=psi[2](0),D(psi[1])(0)=1,psi[1](1)=0};
res:=dsolve(sys); #Analytical solution actually possible
res2:=evalf(res);
plot(subs(res2,[psi[1](x),psi[2](x)]),x=0..1);
resnum:=dsolve(sys,numeric); #Numerical solution
plots:-odeplot(resnum,[[x,psi[1](x)],[x,psi[2](x)]],0..1);

You say that the input should be the differential equation, but you use f inside as a function.
With f as input your procedure should work basically correct, but here is a cleaner version

eul:=proc(f,h,x0,y0,xn)
  local no_points,x1,y1,i;
  no_points:=round(evalf((xn-x0)/h));
  x1:=x0;
  y1:=y0;
 
  for i from 1 to no_points do
      y1:=y1+evalf(h*f(x1,y1));
      x1:=x1+h
  end do;
  y1
end proc;

f:=(x,y)->x^2*y^3; #Input the isolated right hand side as a function of two variables.
eul(f,.01,0,1,1); #Stepsize 0.01.

#Comparison to analytical solution:
res:=dsolve({diff(y(x),x)=f(x,y(x)),y(0)=1});
evalf(eval(rhs(res),x=1));

Here is an illustration of what I believe you are asked to do:

restart;
A:=LinearAlgebra:-RandomMatrix(5,generator=-0.1..0.1);
Student:-NumericalAnalysis:-SpectralRadius(A);
f0:=<1,2,0,5,-6>;
f:=f0:
for i to 15 do f:=A.f end do:
f;

To prove anything you need to know what the spectral radius is. You can consult the help page for
Student:-NumericalAnalysis:-SpectralRadius,
where it is clearly stated.
You may also want to consult
http://en.wikipedia.org/wiki/Spectral_radius
Secondly you need to establish an inequality from

f[n+1]-f[n] = A.( f[n] - f[n-1] ) of this sort:

| f[n+1]-f[n] | <= k* | f[n]-f[n-1] |  (where| | is a norm, and k < 1.)
(In full generality this may have to be modified using Gelfand's formula shown in the Wikipedia reference).
Now applying this repeatedly you get

| f[n+1]-f[n] | <= k^n* | f[1]-f[0] |.

Thus the infinite series Sum( | f[n+1]-f[n] |, n=0..infinity) is convergent. Consequently the series
Sum(  f[n+1]-f[n] , n=0..infinity) is (pointwise) convergent. But it is a telescoping sum so the sequence {f[n]} is convergent (pointwise).

Telescoping means that for all N>=0:   Sum(f([n+1]-f[n],n=0..N) = f[N+1] - f[0]

Is there more to it than

r:=3: L:=12:
plot3d([r*cos(theta),r*sin(theta),z],theta=0..2*Pi,z=0..L);

If you replace the deprecated array with Array then you see something.
Alternatively (if you insist on using array) end the procedure with eval(M).
No need for return in either case, but it doesn't hurt.
If you actually want a Matrix then replace array with Matrix.

Edit: Boldface added for clarity. -Carl Love

You can write the two odes as one.

restart; # Notice that Restart (capital R) has no effect (to catch that use semicolon, not colon)
a:=0.13:
b:=0.41:
reynolds:=1.125*10^8:
Eq1:=diff(f(x),x$3)+diff(f(x),x$2)*f(x)+b^2*sqrt(2*reynolds)*diff(diff(f(x),x$2)*f(x)*x^2,x$1);
Eq2:=diff(g(x),x$3)+diff(g(x),x$2)*g(x)+c*a^2*sqrt(2*reynolds)*diff(diff(g(x),x$2)*x,x$1);
eq1:=isolate(Eq1,diff(f(x),x,x,x));
eq2:=subs(g=f,isolate(Eq2,diff(g(x),x,x,x)));
EQ:=diff(f(x),x,x,x)=piecewise(x<c*0.1,rhs(eq1),rhs(eq2));
#You are basically free to set c as what you want and f '' (0 to what you want.
#So you could just do that. Otherwise use the parameters option:
res:=dsolve({EQ,f(0)=0,D(f)(0)=0,(D@@2)(f)(0)=f2},numeric,parameters=[f2,c],output=listprocedure);
res(parameters=[1,10]); #Setting values for f '' (0) and c (in that order)
plots:-odeplot(res,[seq([x,diff(f(x),[x$i])],i=0..2)],0..2); #This plots from and past 0.1*c

The error message tells the story.

I tried adding F(0)=0, as one condition is missing on F.

bcs2 := f(0) = 0, (D(f))(0) = 1, (D(f))(N) = 0, F(N) = f(N), theta(0) = 1, theta(N) = 0, theta1(N) = 0,F(0)=0;
L := [0, 5, 10];

I had success for M =0 and M=5:
for k to 3 do R[k] := dsolve(eval({Eq1, Eq2, Eq3, Eq4, bcs2}, M = L[k]), [f(eta), F(eta), theta(eta), theta1(eta)], numeric, output = listprocedure,method=bvp[midrich]) end do;

p1:=odeplot(R[1],[[eta,f(eta)],[eta,F(eta)]],0..N):
p2:=odeplot(R[1],[[eta,theta(eta)],[eta,theta1(eta)]],0..N):
display(Array([p1,p2]));
#Try the same for R[2].
If you change N to 5 you get success for all three values of M.

Your function has period Pi and is even. On the interval (-Pi/2,Pi/2) you can remove abs.
Thus f has a cosine expansion of the kind you mention.
f1:=ln(cos(x));
plot(f1,x=-Pi/2..Pi/2);
an1:=4/Pi*Int(f1*cos(2*n*x),x=0..X); #limit x->Pi/2 (left) later
a0:=value(eval(an1,{n=0,X=Pi/2})); #n=0 is special
IntegrationTools:-Parts(an1,f1);
an:=limit(%,X=Pi/2,left) assuming n::posint;
seq(value(eval(an,n=k)),k=1..10); #Concrete values help the integration
S:=a0/2+Sum((-1)^(n+1)/n*cos(2*n*x),n=1..infinity); #Not an unreasonable guess
plot([f1,eval(S,infinity=10)],x=-Pi/2..Pi/2); #Graphical confirmation of the truncated series.
#
#  Now backwards:
S;
#So consider
-ln(2)-Sum((-z)^n/n,n=1..infinity);
value(%);
eval(%,z=exp(2*I*x));
evalc(Re(%));
simplify(%);
expand(%);
simplify(%) assuming cos(x)>0;

I found no problem with either of these:

DE1 := {diff(z(x), x) = y(x)*x, diff(y(x), x, x) = y(x)*z(x)};

DEtools[DEplot3d](DE1, [y(x), z(x)], x = -2 .. 2, [[y(0) = 1, (D(y))(0) = 2, z(0) = 1]], y = -4 .. 4, obsrange = true, stepsize = 0.5e-1, iterations = 3, orientation = [-124, 72]);

DEtools[DEplot](DE1, [y(x), z(x)], x = -2 .. 2, [[y(0) = 1, (D(y))(0) = 2, z(0) = 1]], y = -4 .. 4,scene=[y,z]);

Not nearly as nice as Axel Vogt's result, but an analytic result nevertheless, can be found by using the changes of variables which Maple are using as revealed by setting infolevel[int]:=5.

restart;
A:=Int(cosh(t)/(cosh((17/15)*t)+cosh(t)), t = 0 .. infinity);
B:=simplify(convert(A,exp));
IntegrationTools:-Change(B,t=15/2*ln(x));
IntegrationTools:-Expand(%);
#value(%); #No useful result
r:=int(x^k/(x^16+1),x=1..infinity);
res:=15/2*add((-1)^k*r,k=0..14);
evalf(res);
convert(res,hypergeom);

I hope I got the dynamics right:

restart;
L := -sin(omega*t-phi(t))*(diff(phi(t), t))*a*l*m*omega+(1/2)*(diff(phi(t), t))^2*l^2*m+cos(phi(t))*g*l*m;
L1:=eval(L,{diff(phi(t),t)=phi1,phi(t)=phi});
ode:=diff(eval(diff(L1,phi1),{phi=phi(t),phi1=diff(phi(t),t)}),t)-eval(diff(L1,phi),{phi=phi(t),phi1=diff(phi(t),t)})=0;
params:={g=9.81,omega=1,l=2,a=1,m=1};
eval(ode,params);
res:=dsolve({eval(ode,params),phi(0)=Pi/4,D(phi)(0)=0},numeric,output=listprocedure);
Phi:=subs(res,phi(t));
p:=t->plots:-display(plot([[0,0],[2*sin(Phi(t)),-2*cos(Phi(t))]],color=red,thickness=4),
                     plottools:-disk([2*sin(Phi(t)),-2*cos(Phi(t))],.1,color=blue),
                     scaling=constrained,view=-2.2..0);
p(5);
plots:-animate(p,[t],t=0..10,frames=100);

You can find the laplace transform of the solution.

restart;
pde:=diff(u(x,t),x,x) = diff(u(x,t),t);
res:=inttrans[laplace](pde,t,s);
ode:=eval(res,{u(x,0)=x*(Pi/4-x),laplace(u(x,t),t,s)=w(x)});
L:=dsolve({ode,w(0)=0,w(Pi/4)=0});
convert(rhs(L),trigh);
inttrans[invlaplace](%,s,t);
Maple doesn't find the inverse.