@tomleslie 

Thanks i have seen it, i readlly appreciate. But i have a question, what if i want to define a function to plot tetta' w.r.t episilon

@noah2gud 

Actually, the fsolve solution is nt the same with that of the numerical solution of the same problem.

@nidojan

it was a mistake, i upload a different code.

The real problem i want to solve is thisadm_shateyi.mwadm_shateyi.mw 

@tomleslie 

i use maple 17, bt i could nt get the result for my alpha[1] and alpha[2]

@nidojan 

i need the answer for alpha[1] and alpha[2]

@tomleslie 

thanks. but i still have a problem with the solution, comparing the dtm solution along with the numerical, the answer are not the same.  Any help?

restart;
Digits := 16;
with(ODETools);
with(student);
with(plots);
inf := 9;
equ1 := diff(f(eta), `$`(eta, 3))-a*(diff(theta(eta), eta))*(diff(f(eta), eta, eta))/(1+a*theta(eta))+m(1+a*theta(eta))*(1-(diff(f(eta), eta))^2)+M*(1+a*theta(eta))*(1-(diff(f(eta), eta)))+(1+a*theta(eta))*(diff(f(eta), eta, eta))*f(eta)+(1+a*theta(eta))*lambda*theta(eta) = 0;
equ2 := ((diff(theta(eta), eta))^2*N*b*(1+b*theta(eta))^(N-1)+(diff(theta(eta), eta, eta))*(1+b*theta(eta))^N+4*R/(3*Pr).(diff(theta(eta), eta, eta)))/Pr+Ec*(diff(f(eta), eta, eta))^2-(diff(f(eta), eta))*theta(eta)*(2*m-1)+(diff(theta(eta), eta))*((m+1)*(1/2))*f(eta)+M*Ec*(1-(diff(f(eta), eta)))^2 = 0;
Bcs := f(0) = Fw, (D(f))(0) = `ε`+delta*((D@@2)(f))(0), (D(f))(inf) = 1, (D(theta))(0) = -1, theta(inf) = 0;
NULL;

Pr := 6; M := 1; m := 1; delta := 1; Ec := 1; R := 1; b := 1; Fw := 1; N := 1; lambda := 1; `ε` := 1; a := 0.1e-1;
SolP1 := dsolve({Bcs, equ1, equ2}, numeric);
SolP1(0);
[
[eta = 0., f(eta) = 1.000000000000000,
[

d
----- f(eta) = 1.006562455989112,
deta

d / d \
----- |----- f(eta)| = 0.006562455989111793,
deta \ deta /

theta(eta) = 0.1782538387967585,

d ]
----- theta(eta) = -1.000000000000000]
deta ]

@Carl Love 

restart;
Digits := 5;
with(ODETools);
with(student);
with(plots);
inf := 5;
with(LinearAlgebra);
equ1 := (1+2*n)*f(eta)*(diff(theta(eta), eta))/(1+3*n)-(diff(theta(eta), `$`(eta, 2))) = 0;
equ2 := ((1+n)*(diff(f(eta), eta))^2/(1+3*n)-(1+2*n)*f(eta)*(diff(f(eta), eta, eta))/(1+3*n))/Bo+(diff(f(eta), `$`(eta, 3)))^n-theta(eta) = 0;
Bcs := f(0) = 0, (D(f))(0) = 0, (D(f))(inf) = 0, theta(0) = 1, theta(inf) = 0;
Bo := 1; n := 2;
SolP1 := dsolve({Bcs, equ1, equ2}, numeric);
Error, (in dsolve/numeric/bvp/convertsys) unable to convert to an explicit first-order system
SolP1(0);

@tomleslie 

Thanks for your time, i really appreciate.

when i choose a value for my infinity. i have the value of A and B. but still my result is not the same compare to the numerical solution