Following the line of thought in my answer to your previous question ( http://www.mapleprimes.com/questions/133924-Multiple-Solution-For-Ode) you can do some exploration like this:

restart;
with(plots):
eq:=.7246228659*(diff(f(x), x, x, x))+f(x)*(diff(f(x), x, x))+.6666666666-.6666666666*(diff(f(x), x))^2=0;

b:=10: #Taking infinity to be 10.
#We solve the initial value problem at x = 0 with f''(0) as the parameter p2 and a as another parameter:
res := dsolve({eq = 0, f(0) = 0, D(f)(0) = a, (D@D)(f)(0) = p2}, numeric,output=listprocedure,parameters=[a,p2]);
#The numerical procedures for f, f', and f'':
F0,F1,F2:=op(subs(res,[f(x),diff(f(x),x),diff(f(x),x,x)]));
#Making a procedure p which computes f'(b) - 1, i.e. F1(b) - 1
p:=proc(aa,pp2) if not type([aa,pp2],list(numeric)) then return 'procname(_passed)' end if;
res(parameters=[aa,pp2]);
F1(b)-1
end proc;
#Testing p:
p(-1.2,2);
#Plotting p(-1.1,pp2):
plot(p(-1.1,pp2),pp2=-10..10);
#In the following animation the bounded increasing function tanh is applied to keep values from getting enormous.
animate(plot,[tanh(p(aa,pp2)),pp2=0..1.2],aa=-1.2..-1,frames=25);
#Attempt to narrow down a critical a-value:
for k from 0 to 5 do A:=-1.099-k*.0001; fsolve(p(A,pp2)=0,pp2=0..1.2) end do;
#A 3d-plot:
plot3d(tanh(p(aa,pp2)),aa=-1.2..-1,pp2=0..1.2,axes=boxed);

I have made a couple of corrections: Multiplication signs after f1 and f2 are removed, and pointplot is used.

In order to include the initial point in the plot I have inserted a line and changed the line defining p[i] to one defining p[i+1].

EulerSystem:=proc (f1,f2,t0,x0,v0,h,n) local t,x,v,i,p;
t[0]:=t0; x[0]:=x0; v[0]:=v0;
p[0]:=plots[pointplot]([t[0],x[0]]);
for i from 0 to (n-1) do;
x[i+1]:=x[i]+h*f1(t[i],x[i],v[i]);
v[i+1]:=v[i]+h*f2(t[i],x[i],v[i]);
t[i+1]:=t[i]+h;
p[i+1]:=plots[pointplot]([t[i+1],x[i+1]]);
end do:
plots[display](seq(p[i],i=0..n));
end proc:

 

f1:=(t,x,v) -> v:

f2:=(t,x,v) -> -(0.006549/0.7038)*v-sin(x/1.03)*9.82:

EulerSystem(f1,f2,0,55.1,0,0.05,1000);

#A different approach where the plotting is done outside the procedure:

restart;
EulerSystem:=proc (f1,f2,t0,x0,v0,h,n) local t,x,v,i;
t:=Vector(n+1);
x:=Vector(n+1);
v:=Vector(n+1);
t[1]:=t0; x[1]:=x0; v[1]:=v0;
for i from 1 to n do;
x[i+1]:=x[i]+h*f1(t[i],x[i],v[i]);
v[i+1]:=v[i]+h*f2(t[i],x[i],v[i]);
t[i+1]:=t[i]+h;
end do;
t,x,v;
end proc:

 

f1:=(t,x,v) -> v:

f2:=(t,x,v) -> -(0.006549/0.7038)*v-sin(x/1.03)*9.82:

T,X,V:=EulerSystem(f1,f2,0,55.1,0,0.05,1000);
plot(T,X);

The initial value problem at 0 has a unique solution, your boundary value problem does not.

You could do like this, where I begin by giving Maple the boundary value problem. Since you didn't give us b I chose b = 1.

restart;
with(plots):
eq1 := .6740221708*diff(f(x), x, x, x)+f(x)*diff(f(x), x, x)+1-diff(f(x), x)^2=0;
b:=1:
sol := dsolve({eq1 = 0, f(0) = 0, D(f)(0) = -1.2, D(f)(b) = 1}, numeric,output=listprocedure);
odeplot(sol,[x,f(x)],0..1);
#Saving the plot for later:
PL1:=%:
#The numerical procedures for f, f', and f'':
f0,f1,f2:=op(subs(sol,[f(x),diff(f(x),x),diff(f(x),x,x)]));
f2(0);
#Now we solve the initial value problem at x = 0 with f''(0) as the parameter p2:
res := dsolve({eq1 = 0, f(0) = 0, D(f)(0) = -1.2, (D@D)(f)(0) = p2}, numeric,output=listprocedure,parameters=[p2]);
#The numerical procedures for f, f', and f'':
F0,F1,F2:=op(subs(res,[f(x),diff(f(x),x),diff(f(x),x,x)]));
#Testing with p2 = f2(0) found initially:
res(parameters=[f2(0)]);
res(b);
odeplot(res,[x,f(x)],0..1);
F1(b);
#Making a procedure p which computes f'(b) - 1, i.e. F1(b) - 1
p:=proc(pp2) if not type(pp2,numeric) then return 'procname(pp2)' end if;
res(parameters=[pp2]);
F1(b)-1
end proc;
#Testing p:
p(f2(0));
#Plotting p:
plot(p,-20..10);
#In the interval [-20,10] there are apparently 2 zeros. The one we knew and a new one near -10:
pm2:=fsolve(p(pp2),pp2=-10);
res(parameters=[pm2]);
F1(b);
#The new solution:
odeplot(res,[x,f(x)],0..1,color=blue);
PL2:=%:
#Both solutions
display(PL1,PL2);

Try wrapping the definitions of f1 and g1 in 'evalc' and omit the x- and y- ranges in DEplot:

restart;
with(DEtools):

z:=exp(I*((t/(2*Pi))+(2*Pi*j/3)));

f1:=evalc(Re((I/(2*Pi))*((1/(subs(j=1,z)-(x(t)+I*y(t))))+(1/(subs(j=2,z)-(x(t)+I*y(t))))+(1/(subs(j=3,z)-(x(t)+I*y(t)))))));

g1:=evalc(-1*Im((I/(2*Pi))*((1/(subs(j=1,z)-(x(t)+I*y(t))))+(1/(subs(j=2,z)-(x(t)+I*y(t))))+(1/(subs(j=3,z)-(x(t)+I*y(t)))))));

sys1:=[diff(x(t),t)=f1,diff(y(t),t)=g1]:

DEplot(sys1, [x(t), y(t)], t = -20 .. 10 ,  [[x(0) =0, y(0) = 0]],linecolor = blue, thickness = 2, stepsize = .1, arrows = medium, color = red);

Change the variable of integration from phi to e.g. phi1. The variable of integration is not local to 'int', but is in your definition considered the formal variable used in defining the procedure. (And use 'Pi', not 'pi').

E := phi->int(1/(1-(1-1/kappa^2)*sin(phi1)^2)^(1/2), phi1 = 0 .. phi);

E(Pi/2);

The solution reaches infinity in finite time as can also be seen from the exact solution, which can be obtained in this case:

dsolve({eq2, x(0) = 1, y(0) = -1});

Y:=subs(%,y(t));
fsolve(1/Y=0,t=1.3);
fsolve(1/Y=0,t=-1.1);
plot(Y,t=-2..2,-10..10);

You could differentiate the equations w.r.t. t as in the following example done in Maple 12. If you start with 'restart' then p1 and p2 are reproducible.

restart;
p1:=randpoly([x,y,dx,dy,t]);
p2:=randpoly([x,y,dx,dy,t]);
eq1:=subs(x=x(t),y=y(t),dx=diff(x(t),t),dy=diff(y(t),t),p1);
eq2:=subs(x=x(t),y=y(t),dx=diff(x(t),t),dy=diff(y(t),t),p2);
#Notice that eq1 really means eq1 = 0.
#This fails:
dsolve({eq1,eq2,x(0)=1,y(0)=2},numeric);
deq1:=diff(eq1,t);
deq2:=diff(eq2,t);
#Finding initial values for the derivatives:
eval({p1,p2},{t=0,x=1,y=2});
dxdy:=fsolve(%,{dx,dy});
#There may be other solutions, so here is a job that remains to be done!
#Now solving the second order system:
res:=dsolve({deq1,deq2,x(0)=1,y(0)=2,D(x)(0)=subs(dxdy,dx)  ,D(y)(0)=subs(dxdy,dy) },numeric);
plots:-odeplot(res,[[t,x(t)],[t,y(t)]],0..1);

Besides real problems there are the following errors:

'pi' should be 'Pi'.

The fifth derivative of the a's should be excluded from the initial conditions.

The syntax in the definition of 'dsolved' should be

dsolve(SystemOfEq union inits,[a_1(t),a_2(t),a_3(t),Omega(t),Lambda(t)],numeric);

(or SytemOfEq as you spell it)

with additional options like 'stiff=true'.

It is not quite clear to me what you want, but try this:

restart;
m[1] := 30; c[1] := 48000; c[0] := 20000; d[1] := 21; d[0] := 70;
#S, DS, and DDS are taken from your previously posted question:
S:=9.28+14.97*sin(t)-14.29*cos(t)-.11*sin(15*t)-0.2e-1*cos(15*t)+0.4e-1*sin(14*t)-0.5e-1*cos(14*t)+0.4e-1*sin(13*t)-0.7e-1*cos(13*t)+.11*sin(12*t)+.16*cos(12*t)-.1*sin(11*t)+0.4e-1*cos(11*t)-.27*sin(10*t)+0.7e-1*cos(10*t)+.13*sin(9*t)-.61*cos(9*t)+.52*sin(8*t)+.28*cos(8*t)+.13*sin(7*t)+.78*cos(7*t)-1.76*sin(6*t)-.32*cos(6*t)+.88*sin(5*t)-1.21*cos(5*t)+.47*sin(4*t)-.26*cos(4*t)+1.29*sin(3*t)+3.16*cos(3*t)-12.12*sin(2*t)-.89*cos(2*t);
DS:=diff(S,t);
DDS:=diff(DS,t);
DGL1 := m[1]*(D(D(x)))(t)+d[0]*(D(x))(t)+d[1]*(D(x))(t)+c[0]*x(t)+c[1]*x(t)+50000 = -(1/1000)*c[1]*S-(1/1000)*d[1]*DS*((1/60)*(2*Pi*250))-(1/1000)*m[1]*DDS*((1/60)*(2*Pi*250))^2;
init := x(1.1) = 1, (D(x))(1.1) = .2;
A := dsolve({DGL1, init},x(t),method=laplace):
#Plotting the solution:
plot(rhs(A),t=0..5);
# and on a smaller interval:
plot(rhs(A),t=1.1..1.2);
fsolve(rhs(A)=0,t=1.12);

#When you say that it is required that the left hand side of the equation be positive, it is not clear to me which equation you are talking about, in view of your previously posted question.

You could use DEplot in the DEtools package, or you could use dsolve/numeric and then plot the result using odeplot from the plots package.

It may be the same problem as the one you get when you try the following:

X:=x->sin(x);
Y:=x->ln(x);
#X and Y are procedures, but arctan doesn't accept a procedure as input
plot(arctan(X/Y),2..5);
#Same problem, but more directly:
plot(arctan(sin/ln),2..5);
#In this case just do:
plot(arctan(X(s)/Y(s)),s=2..5);

I found it interesting that it makes a difference if you exclude the tautology D(y)(0) = D(y)(0) from the conditions:

restart;
eq:=diff(y(x),x$4)+(2/x)*diff(y(x),x$3)-((2*(n-1)^2+1)/(x^2))*diff(y(x),x$2)+((2*(n-1)^2+1)/(x^3))*diff(y(x),x)+(((n-1)^4-4*(n-1)^2)/(x^4))*y(x)=8*x^(n-1)*(n-n*x^2-x^2)+(n+1)*(n+3)*(n-3)*(n-5)/12-(n+2)*(n+4)*(n-4)*(n-6)*x/12;
res1:=dsolve({eq,y(1)=0,D(y)(1)=0,D(y)(0)=D(y)(0)});
res2:=dsolve({eq,y(1)=0,D(y)(1)=0});
indets(res1,name);
indets(res2,name);

#To proceed from res1 you could do as follows:

combine(rhs(res1),power);
#Just inspect the output and copy the following requirement:
24*n^2*_C5-n^2+72*n*_C5+3*n+48*_C5+16=0;
# and proceed as follows:
solve(%,{_C5});
res:=eval(res1,%);

#It is reassuring that if you proceed from res2 you obtain the same result:

combine(rhs(res2),power);
#Now there are two requirements:
24*_C3*n^3+24*_C3*n^2+24*_C4*n^3+3*n^2-96*n*_C3+48*_C4*n^2+5*n-96*_C3-24*n*_C4-4-48*_C4=0;
-n*_C3+_C3-n*_C4-1/8=0;
solve({%,%%},{_C3,_C4});
resN:=eval(res2,%);
rhs(res)-rhs(resN);
simplify(%);

Here is a very fast and straightforward approach:

restart;
with(plots) :
mm := 3: nn := 3: a := 5.: e := 1.:  R := .1: h := 1.:

for n to nn do
for nd to nn do
for md to mm-1 do
integ[n,nd,md]:=evalf(Int( sin(Pi*n*x/a)*sin(Pi*nd*x/a)*sin(Pi*n*y/a)*sin(Pi*nd*y/a)/sqrt((x-y)^2+4*R^2*sin(Pi*md/mm)^2), y = 0. .. a, x = 0. .. a))
end do
end do
end do;

sys:={seq(diff(Y[n](t),t)=-I*e^2*4/(h*a^2)*Y[n](t)*add(add(abs(Y[nd](t))^2*integ[n,nd,md], md = 1 .. mm-1),   nd = 1 .. nn),n=1..nn)};
inits:={seq(Y[n](0)=1/sqrt(nn), n = 1 .. nn)};
res:=dsolve(sys union inits,numeric,range=0..0.1,method=rkf45);
odeplot(res,[seq([t,abs(Y[n](t))],n=1 .. nn)]);

You may try this (done in Maple 12, as it so happens).

sol := dsolve(Systeme, type = numeric,range=0..20,output=listprocedure);

 odeplot(sol, [T1(t), T2(t)], 0 .. 20,refine=1);

TT1,TT2:=op(subs(sol,[T1(t),T2(t)]));
p:=t->pointplot([cos(TT1(t)), sin(TT1(t))],color=red,symbol = solidcircle, symbolsize = 25):
animate(p,[t],t=0..8,frames=300);

# Added later: In Maple 15 (as I suspected) you can just do

animate(pointplot,[[cos(TT1(t)), sin(TT1(t))],color=red,symbol = solidcircle, symbolsize = 25],t=0..8,frames=300);

#The need for defining the procedure 'p' above is due to a weakness in 'animate' present in Maple versions earlier than 15, but corrected in Maple 15.

Replace Epr11 := evalc(conjugate(psi1*psi1));

by

Epr11 := evalc(  conjugate(psi1(x))   *  psi1(x)  );

and similarly for the other quantities.