@Amare dsolve/numeric doesn't return a module.
If sol2:=dsolve(....,numeric); then by default a procedure is returned (a procedurelist).
If sol2:=dsolve(....,numeric,output=listprocedure); then a list of procedures is returned.

@Amare dsolve/numeric doesn't return a module.
If sol2:=dsolve(....,numeric); then by default a procedure is returned (a procedurelist).
If sol2:=dsolve(....,numeric,output=listprocedure); then a list of procedures is returned.

@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!).