@Rouben Rostamian  Before trying to solve the problem of this topic, I had tried a simply-supported beam, with successful results. Regarding the problem you suggest, as you previewed, Maple cannot handle it, too.

Worksheet2.mw

@Carl Love When I apply this comand, Maple actually creates two independent plots. I can export them indivudually, but, for some reason, when I try to export the colorbar in EPS, Maple does not export it right (both images are shown below). I hope that, if I get to join both plots, Maple might export it correctely. Alternatively, I would be satisfied with correctly exporting the colorbar, so that I can assemble both graphics in another software.

Thank you for all the answers; both solutions were helpful. In case if helps other people, the code which was enough for my problem was:

restart:

omega0:= 10/3:
plot(omega^2/omega0, omega = 1..1.5,
legend = [typeset(`ω`[0] = evalf[4](omega0), " rad/s")]);

@Rouben Rostamian Thank you very much! I'll keep the work based on the algorythm you have made. 

@Rouben Rostamian  

Thanks for your answer.

1) I wanted to start with a simple ODE because this way I could verify the confiability of my algorythm. The equation I'm really interested in solving is:

x'' + A*x + B*x*(x')^2 + C*x^3 = 0, where A,B and C are known parameters.

The scaling parameter ("e") is applied to the initial conditions, which should have the form:

x(0) = e*x[1](0) + e^2*x[2](0) + h.o.t

x'(0) = diff(x,t) | x = 0

2) You're right. x[j] depends on T[0],T[1] and T[2]; not only on T[1].

3) I agree with you, but, since I will need to use the MMS for a branch of equations, I thought the automatising would be a good idea. I have done it manually for a first approximation. Doing a second one, though, with terms of higer order of "e", would be easier if done autommatically. The main issue here is to rewritte the ODE, using the variable transformations suggested by the MMS.