You made it very clear but what will be the possible staring point to handle the situation?

@PatrickT, I tried to follow ur suggestions but no luck.

1. complex=false

 dsolve(subs(para,epsilon=1.6, {bc, eq1, eq2}), numeric,complex=false):
Error, (in dsolve/numeric/bvp) bvp keyword was complex, optional keyword must be

 one of 'output', 'value', 'range', 'abserr', 'maxmesh', 'initmesh','continuation', 'mincont', 

'approxsoln', 'adaptive', 'interpolant', 'optimize'

 2. Re

eq1:=Re(n*(-diff(f(eta),eta$2))^(n-1)*diff(f(eta),eta$3))+(2*n/(n+1))*f(eta)*

diff(f(eta),eta$2)-diff(f(eta),eta$1)^2-M*diff(f(eta),eta$1)+M*epsilon+epsilon^2+

lambda*theta(eta)=0;

eq2:=remain the same with bcs;

dsolve(subs(para,epsilon=1.6, {bc, eq1, eq2}), numeric):
Error, (in fproc) unable to store 'HFloat(0.251188643150958)*

RootOf(2000000000*Re(_Z)-HFloat(6.5e10))' when datatype=float[8]

Thanks

@PatrickT, I tried to follow ur suggestions but no luck.

1. complex=false

 dsolve(subs(para,epsilon=1.6, {bc, eq1, eq2}), numeric,complex=false):
Error, (in dsolve/numeric/bvp) bvp keyword was complex, optional keyword must be

 one of 'output', 'value', 'range', 'abserr', 'maxmesh', 'initmesh','continuation', 'mincont', 

'approxsoln', 'adaptive', 'interpolant', 'optimize'

 2. Re

eq1:=Re(n*(-diff(f(eta),eta$2))^(n-1)*diff(f(eta),eta$3))+(2*n/(n+1))*f(eta)*

diff(f(eta),eta$2)-diff(f(eta),eta$1)^2-M*diff(f(eta),eta$1)+M*epsilon+epsilon^2+

lambda*theta(eta)=0;

eq2:=remain the same with bcs;

dsolve(subs(para,epsilon=1.6, {bc, eq1, eq2}), numeric):
Error, (in fproc) unable to store 'HFloat(0.251188643150958)*

RootOf(2000000000*Re(_Z)-HFloat(6.5e10))' when datatype=float[8]

Thanks

@Preben Alsholm thanks for your explanation.

The system is exactly the one for which I need a solution. You can see the attached research

paper Magnetohydrodynamic_.pdf, where they give the solution of exactly the same equation (eq1).

Infact, I adopted the eq1 from there including one extra term which I needed for my model. 

As you rightly mentioned that,  "if diff(f(eta),eta,eta) is real and positive then (-diff(f(eta),eta$2))^(n-1) is imaginary". What I am thinking is that can we consider the real part and find the solution?

How we will do it in Maple?

Thanks 

@Preben Alsholm thanks for your explanation.

The system is exactly the one for which I need a solution. You can see the attached research

paper Magnetohydrodynamic_.pdf, where they give the solution of exactly the same equation (eq1).

Infact, I adopted the eq1 from there including one extra term which I needed for my model. 

As you rightly mentioned that,  "if diff(f(eta),eta,eta) is real and positive then (-diff(f(eta),eta$2))^(n-1) is imaginary". What I am thinking is that can we consider the real part and find the solution?

How we will do it in Maple?

Thanks 

I am not sure about it but here is the thing

Its an initial condition of that particular type "D(f)(0)=beta/(1-phi)^(2.5)*(D@D)(f)(0)".

When I try to solve the ode I do not need to specify (D@D)(f)(0) to get the

desired results. 

bc:=f(0)=0,D(f)(0)=beta/(1-phi)^(2.5)*(D@D)(f)(0),D(f)(N)=1;

N:=5:beta:=1:M:=1:phi:=0.1:

d:=dsolve({eq1,bc},numeric,output=array([seq( 0.01*i, i=0..100*N)])):

d1:=odeplot(d, [[eta,diff(f(eta),eta), color=red]],0..N):

display(d1);

One thing more M must be +ve, the same is for beta, M belongs to [0, 3].

I am not sure about it but here is the thing

Its an initial condition of that particular type "D(f)(0)=beta/(1-phi)^(2.5)*(D@D)(f)(0)".

When I try to solve the ode I do not need to specify (D@D)(f)(0) to get the

desired results. 

bc:=f(0)=0,D(f)(0)=beta/(1-phi)^(2.5)*(D@D)(f)(0),D(f)(N)=1;

N:=5:beta:=1:M:=1:phi:=0.1:

d:=dsolve({eq1,bc},numeric,output=array([seq( 0.01*i, i=0..100*N)])):

d1:=odeplot(d, [[eta,diff(f(eta),eta), color=red]],0..N):

display(d1);

One thing more M must be +ve, the same is for beta, M belongs to [0, 3].

@Joe Riel

One thing more, I am trying to get the data in .txt file for the four plots using this inside the loop

MM[M]:=

;

ExportMatrix("C:/data.txt", MM[M]):

But i just get data for one plot. How to do it for four of them.

@Joe Riel

One thing more, I am trying to get the data in .txt file for the four plots using this inside the loop

MM[M]:=

;

ExportMatrix("C:/data.txt", MM[M]):

But i just get data for one plot. How to do it for four of them.

Thanks @Joe Riel, it worked. But how to use different colors for the different output.

Thanks @Joe Riel, it worked. But how to use different colors for the different output.

Well, thanks guys (@Preben Alsholm & @herclau) for your helpful suggestions.

Actually, I am looking for some thing/example which clearly show me

the change in %. Some numerical/graphical thing/example, where we can see the

change.

Thanks

Well, thanks guys (@Preben Alsholm & @herclau) for your helpful suggestions.

Actually, I am looking for some thing/example which clearly show me

the change in %. Some numerical/graphical thing/example, where we can see the

change.

Thanks

It seems that Maple's numerical method for solving boundary-value problems

can only return one solution at a time. Therefore, when lambda < roughly -1,

it can only give you the value of -theta'(0) corresponding to the "lower branch solution" ,

because that is the solution to which the numerical solver has converged.

I think most simple methods are iterative; it may then be possible to give it an initial guess

for the solution and by choosing that guess carefully to converge on the upper rather than

the lower-branch solution.

For this purpose I have simplified the problem under consideration and converted it into IVP

restart:
with(plots):

eq1:= diff(f(eta), eta, eta, eta)+f(eta)*(diff(f(eta), eta, eta))-(diff(f(eta), eta))^2

+1+k*(f(eta)^2*(diff(f(eta), eta, eta, eta))-2*f(eta)*(diff(f(eta), eta))*(diff(f(eta), 

eta, eta))) = 0;

eq2:=diff(theta(eta), eta$2)+pr*f(eta)*diff(theta(eta), eta$1)=0;

respar:=dsolve({eq1,eq2,f(0)=s,D(f)(0)=lambda,theta(0)=1,(D@D)(f)(0)=f2/(1+k),

D(theta)(0)=-f3},numeric,output=listprocedure,parameters=[k,pr,s,lambda,f2,f3]);

X1,X2,F,F1,F2:=op(subs(respar,[theta(eta),diff(theta(eta),eta),f(eta),diff(f(eta),eta),

diff(f(eta),eta,eta)]));

p := proc (N,k,pr,s,lambda,f2,f3) if not type([args], list(numeric)) then return

'procname(args)' end if; respar(parameters = [k,pr,s,lambda,f2,f3]); X1(N): F1(N)-1

end proc;

implicitplot(p(4, 0.5, 0.7, 0.5,lambda, 3.142141596,f3) =0, lambda=-1.75..0,

f3 = 0..7, grid = [35, 35], gridrefine = 2, resolution = 1000, crossingrefine = 3,

color = black);

But as you can see, the output does not make any sense.

Thanks

The idea is to plot the expression in F3 vs x=0..1.

For y we can use output = array([seq(0.1e-1*i, i = -110 .. 110)]).

But I can not make it happen.

Thanks