In the Allan Wittkopf homepage you link to there is a page describing the Rif package. There it says,

"It is the successor of the Standard Form project, and has the following additional features: ...."

My guess is that the commands and features there are available in DEtools.

In Maple 18 try
?Rif

ODE:= diff(y(x),x$2) - 1.0325*diff(y(x),x) + 1.36*y(x) = sin(2*x):
bcs:={y(0)=0, y(1)=1};
#Now introduce the equation suggested by acer in a comment above:
sys:={ODE,diff(y(x),x,x)=w(x)};
T:=DEtools[rifsimp](sys);
T[Solved];
res:=dsolve({T[Solved][]} union bcs,numeric);
plots:-odeplot(res,[x,w(x)]);

@rwether
Hints:
allvalues works on a RootOf expression and returns the roots when possible, which it clearly is here.
evalc(expr) returns expr in the form a+I*b. If there are any names (variables) in expr, evalc assumes that they represent real numbers. 
evalf turns real or complex numbers into floats (decimal numbers).

The order in which I used these followed my thinking: I need the four roots, so allvalues. I see that exponential functions with imaginary arguments are involved. For a problem clearly having a real solution, I would like to see it as such, therefore evalc.
Finally, since the roots are kind of ugly, I use evalf to get decimal numbers.

As it turns out, if you don't want to see the intermediate steps, you can just do
resRootOf:=eval(res,{data});
evalc(evalf( resRootOf ));

That this is sufficient must be because the procedure `evalf/Sum` calls some other procedure (allvalues?), which evaluates the roots of the quartic.

To see that something like that is going on you could try executing these lines:
RootOf(_Z^4+8*_Z^2+5); #returns unevaluated
evalf(RootOf(_Z^4+8*_Z^2+5)); #returns only one of the roots
allvalues(RootOf(_Z^4+8*_Z^2+5)); #Now we have all of the roots
evalf(allvalues(RootOf(_Z^4+8*_Z^2+5)));
#########
Functions.
If you really want functions you can use unapply on each element of the list Z, either one at a time:
f1:=unapply(Z[1],t); #The first, i.e. z1
or all in one line using elementwise application of unapply:
Zfus:=unapply~(Z,t);
Now Zfus[1] is actually exactly the same function as f1.
I am wondering though, why you want the result as functions and not expressions in t.

Actually I would have used the simple way mentioned by Rouben Rostamian above.
With the present example there doesn't seem to be any reason not to go with that, or does there?
Here is a simple test.

restart;
ODE:= diff(y(x),x$2) - 1.0325*diff(y(x),x) + 1.36*y(x) = sin(2*x):
bcs:={y(0)=0, y(1)=1};
SE:=dsolve({ODE} union bcs); #The exact solution.
S1:=dsolve({ODE} union bcs,numeric);
#Differentiating the ODE and finding an additional condition:
ODE3:=diff(ODE,x);
bcs0:=subs(x=0,bcs,convert(ODE,D));
S:=dsolve({ODE3,bcs0} union bcs,numeric);
Y2:=solve(ODE,diff(y(x),x,x));
plots:-odeplot(S1,[x,Y2-diff(rhs(SE),x,x)]);
plots:-odeplot(S,[x,Y2-diff(rhs(SE),x,x)]);

The errors plotted are similar in magnitude.

One might also ask how well is ODE satisfied by S?
plots:-odeplot(S,[x,(lhs-rhs)(ODE)]); #Not bad at all
Clearly the similar question for S1 is not meaningful.

@Met28 You cannot solve the problem via Maple. It is mathematics. The improper integral is simply not convergent. I illustrated that by my Maple statements.

Using the usual definition of convergence in Calculus, you need to show that the improper integrals over the 3 intervals
0..10.69,  10.69..b, b..infinity (where b is chosen so that 10.69 < b < infinity)
are all 3 convergent.

Thus to prove divergence you only need to show that one of the 3 is divergent. As it happens, they are divergent all 3.

Divergence of the first can be proven like this:
res:=int((x*0.187549975)/(x^2-10.69^2), x=0..a) assuming a>0,a<10.69;
limit(res,a=10.69,left);

@Wackeraf I executed all the statements and then examined the odes only. I didn't check the physics, i.e.  whether your derivation of the 4 odes La1, La2, La3, La4 is correct or not. I took them as your code gave them.
That you cannot solve for the highest order derivatives can sometimes be solved by differentiating one of the equations, i.e. using non-algebraic means. It is not obvious to me that that is a way out here.

###
Notice, that in A actually two of the rows are the same:

A[1,..];
A[4,..];

Did you try computing the integral before attempting minimization?

int(0.1e-3+.5*t+0.2e-2*t^2-b*t-a, t = 0 .. 300);

Output:

@acer I tried changing initial condition in your version as I did in my comment above. It doesn't work.

Just checking existing conditions doesn't work either:

ONE_SOLUTION(initial);

The initial condition for D(h) can be included in the parameters, however:

ODE_SOLUTION := dsolve({ODE, h(0) = 0,D(h)(0)=y1}, numeric, method=rosenbrock_dae,parameters = [y1,A, B, C, E, F]):
ODE_SOLUTION('parameters'=[Pi/2,0,0,1,1,1]):

@ There is ambiguity in determining the initial condition even in the case you consider.
The initial condition for u, i.e. u(0)=u1, satisfies cos(u1)=0. Thus there is a choice to be made.
After having found the solution as given, we can do

sol(initial=[0,0,-Pi/2]);

and when plotting we find that we have another solution.

@Carl Love I couldn't agree more!

@Alejandro Jakubi Thank you very much, I shall try to contact Gaston Gonnet.

@Markiyan Hirnyk Yes, but you are using two calls to eval.
You cannot do
eval(f(a*x), x = a, a = 9);

@acer Using method=nonlinearsimplex worked in Maple 15 on my simple example, which is of the same kind as the one described by joha.

@acer I noticed that joha uses Maple 14. On a simple example I tested, this device doesn't work in Maple 15 (and thus most likely not in Maple 14). It does work in Maple 18 (and in Maple 16).

I think it might be necessary to use the operator form with a user supplied gradient as shown by you in 2010 (July 14).

@Benrath plot3d(x*y,x=y..1,y=0..1);