@matteoziff To me the most natural thing to do would be to rewrite the ode as a first order system (sys below) in the standard way using the new additional dependent variable y given by y(t) = diff(x(t),t). Then after that use DEplot from the DEtools package to make a direction field

DEtools[DEplot](sys,[x(t),y(t)], t=0..1, x = -2.1..2.1, y = -1500..1500);

and show that together with the plot from odeplot using display.

But it seems that you don't want to do that?

@Noreen cute I tried the worksheet in Maple 12 (I don't have a working version of Maple 13).

The problem is with fsolve. The syntax used works in Maple 16, but not in Maple 12. I changed that.

So here is a revised version, which works in Maple 12 and I'm sure also in Maple 16.

MaplePrimes12-08-23.mw

@Noreen cute I tried the worksheet in Maple 12 (I don't have a working version of Maple 13).

The problem is with fsolve. The syntax used works in Maple 16, but not in Maple 12. I changed that.

So here is a revised version, which works in Maple 12 and I'm sure also in Maple 16.

MaplePrimes12-08-23.mw

Since this problem is quite similar to the one you posed yesterday:

http://www.mapleprimes.com/questions/136745-How-To-Calculate-Dual-Solution-Problem

I tried the same method. It works fine with b = 6, but certainly failed for b = 20.

But the direct approach you are using also works with b = 6 (even b=6.5). When Pr = 2 you may need something like
sol:=dsolve(eval({de1},Pr=2),numeric,continuation=k);

which worked for me. It starts solving the problem with k=0 and then works its way up to k=1.

@Noreen cute What kind of error?

You could try to execute the following worksheet in Maple 16:

MaplePrimes12-08-23S.mw

@Noreen cute What kind of error?

You could try to execute the following worksheet in Maple 16:

MaplePrimes12-08-23S.mw

Could you tell us how you defined q[n] and how you subsequently defined p[n]?

@samiyare When using the format Spline(X,Y,x) X and Y are not expected to be matrices. So you can just do
X:=convert(X,Vector);
Y:=convert(Y,Vector);
and then continue after that.

@samiyare When using the format Spline(X,Y,x) X and Y are not expected to be matrices. So you can just do
X:=convert(X,Vector);
Y:=convert(Y,Vector);
and then continue after that.

If the result Y found above is close to being a solution, then the solution is not unique. The Y found is not symmetric with respect to x=L/2. But it is immediately seen that if y(x) is a solution to the bvp problem then so is y(L-x). Unless therefore y(L-x) = y(x) for all x in 0..L the problem has more than one solution.

Instead of using fsolve it seems to be a good idea to minimize y(L)^2 + y'(L)^2 using LSSolve from the Optimization package. It is fast and gives great results:

with(Optimization);
res:=LSSolve([p1,p2]);
dsolpar(parameters=convert(res[2],list));
Y(L),Y1(L);
plots:-odeplot(dsolpar,[x,y(x)],0..L);
#Using as initialpoint for LSSolve the result found by fsolve:
res:=LSSolve([p1,p2],initialpoint=[0.7898528584778143356e-2, -0.1372579321779395876e-3]);
dsolpar(parameters=convert(res[2],list));
Y(L),Y1(L);
plots:-odeplot(dsolpar,[x,y(x)],0..L);

If the result Y found above is close to being a solution, then the solution is not unique. The Y found is not symmetric with respect to x=L/2. But it is immediately seen that if y(x) is a solution to the bvp problem then so is y(L-x). Unless therefore y(L-x) = y(x) for all x in 0..L the problem has more than one solution.

Instead of using fsolve it seems to be a good idea to minimize y(L)^2 + y'(L)^2 using LSSolve from the Optimization package. It is fast and gives great results:

with(Optimization);
res:=LSSolve([p1,p2]);
dsolpar(parameters=convert(res[2],list));
Y(L),Y1(L);
plots:-odeplot(dsolpar,[x,y(x)],0..L);
#Using as initialpoint for LSSolve the result found by fsolve:
res:=LSSolve([p1,p2],initialpoint=[0.7898528584778143356e-2, -0.1372579321779395876e-3]);
dsolpar(parameters=convert(res[2],list));
Y(L),Y1(L);
plots:-odeplot(dsolpar,[x,y(x)],0..L);

You could use :-line. By the command with(plottools); the variable line has been defined to be plottools:-line.
Thus you would write
DEplot(deq1,x(t),t=-.5..4,x=-0.2..2,IC,linecolor=black,thickness=3,arrows=:-line);

@JJames _Z1 (printed with a trailing tilde meaning assumptions made on _Z1) stands for an arbitrary integer. So if 'res' above contains _Z1 you can replace it with any integer e.g. 1. 'res' actually also contains a name _C3 (without a trailing tilde). A name like that stands for an arbitrary complex constant.

eval(res,{_Z1=1,_C3=I/2});

As in solving the eigenvalue problem for a matrix A: Find the numbers lambda s.t. Av = lambda*v has a nontrivial solution, there may be ways of finding lambda before actually finding the corresponding v's. In this case we solve the characteristic equation det(A-lambda*Id) = 0.

Above we found the possible values for lambda (actually assuming that they were positive). There were infinitely many, and I chose to use a simpler name n instead of _Z1 or _Z2 or what have you (Maple changes the number if you do it again without a restart).

lambda = (a^2+Pi^2*n^2)^(3/2)/a.

After having found the possible values of lambda you can find all the corresponding eigenfunctions:

dsolve(eval({ode1,ode2,bcs},lambda=(a^2+n^2*Pi^2)^(3/2)/a)) assuming n::integer;

However, Maple can (in principle) solve the whole problem in one command as I noticed in the corrected version of my first answer to your question:

dsolve({ode1,ode2,bcs},{alpha(z),H(z),lambda});

The output, however, is rather unwieldy, but hopefully correct.

With a little more work it can be seen that the corresponding sequence of negative values
lambda= - (a^2+n^2*Pi^2)^(3/2)/a)
are also eigenvalues (n<>0 though!).

@JJames _Z1 (printed with a trailing tilde meaning assumptions made on _Z1) stands for an arbitrary integer. So if 'res' above contains _Z1 you can replace it with any integer e.g. 1. 'res' actually also contains a name _C3 (without a trailing tilde). A name like that stands for an arbitrary complex constant.

eval(res,{_Z1=1,_C3=I/2});

As in solving the eigenvalue problem for a matrix A: Find the numbers lambda s.t. Av = lambda*v has a nontrivial solution, there may be ways of finding lambda before actually finding the corresponding v's. In this case we solve the characteristic equation det(A-lambda*Id) = 0.

Above we found the possible values for lambda (actually assuming that they were positive). There were infinitely many, and I chose to use a simpler name n instead of _Z1 or _Z2 or what have you (Maple changes the number if you do it again without a restart).

lambda = (a^2+Pi^2*n^2)^(3/2)/a.

After having found the possible values of lambda you can find all the corresponding eigenfunctions:

dsolve(eval({ode1,ode2,bcs},lambda=(a^2+n^2*Pi^2)^(3/2)/a)) assuming n::integer;

However, Maple can (in principle) solve the whole problem in one command as I noticed in the corrected version of my first answer to your question:

dsolve({ode1,ode2,bcs},{alpha(z),H(z),lambda});

The output, however, is rather unwieldy, but hopefully correct.

With a little more work it can be seen that the corresponding sequence of negative values
lambda= - (a^2+n^2*Pi^2)^(3/2)/a)
are also eigenvalues (n<>0 though!).

@JohnJames There is no ready made procedure, so yes, you will hve to do it manually. The individual manual steps could, however, be done in Maple, but you will most likely find it much easier to do them on paper.

To illustrate that last point I can give you a few steps, which will probably convince you:

res1:=dsolve(odes[1]) assuming _c[1]<0; #if that is the relevant assumption
res2:=dsolve(odes[2]);
eval(IBC,f=unapply(rhs(op(1,res)),x1,x2));
#We want a nontrivial solution, so
eval(%,{_F1(x1)=1,_F2(x2)=1});
diff(res1,x1);
eval(%,x1=0);
eval(rhs(%%),{x1=L1,_C1=0});
solve(%=0,{_c[1]},AllSolutions);

etc.

In other words, you will have to know what to do yourself at each step. Doing it in Maple just adds to the difficulty in this case.