Question: how to calculate dual solution problem

restart;
b := 7; de1 := diff(f(eta), eta, eta, eta)-M^2*(diff(f(eta), eta))-(diff(f(eta), eta))^2+m1*f(eta)*(diff(f(eta), eta, eta)) = 0, diff(theta(eta), eta, eta)+m1*Pr*f(eta)*(diff(theta(eta), eta)) = 0, f(0) = S, (D(f))(0) = -1, (D(theta))(0) = -1, (D(f))(b) = 0, theta(b) = 0; d1 := subs(m1 = 1, S = 1.75, Pr = 0.7, M = 2, [de1]); d2 := subs(m1 = 1, S = 1.80, Pr = 0.7, M = 2, [de1]); d3 := subs(m1 = 1, S = 1.85, Pr = 0.7, M = 2, [de1]); d4 := subs(m1 = 1, S = 1.90, Pr = 0.7, M = 2, [de1]); da1 := dsolve(d1, numeric); da2 := dsolve(d2, numeric); da3 := dsolve(d3, numeric); da4 := dsolve(d4, numeric); da1(0); da2(0); da3(0); da4(0); with(plots); p1 := odeplot(da1, [eta, diff(f(eta), eta)], linestyle = 1, color = black, legend = "Le = 0.8"); p2 := odeplot(da2, [eta, diff(f(eta), eta)], linestyle = 2, color = black, legend = "Le = 1.0"); p3 := odeplot(da3, [eta, diff(f(eta), eta)], linestyle = 3, color = black, legend = "Le = 1.5"); p4 := odeplot(da4, [eta, diff(f(eta), eta)], linestyle = 4, color = black, legend = "Le = 2.0"); plots[display]({p1, p2, p3, p4}); p5 := odeplot(da1, [eta, theta(eta)], linestyle = 1, color = black, legend = "Le = 0.8"); p6 := odeplot(da2, [eta, theta(eta)], linestyle = 2, color = black, legend = "Le = 1.0"); p7 := odeplot(da3, [eta, theta(eta)], linestyle = 3, color = black, legend = "Le = 1.5"); p8 := odeplot(da4, [eta, theta(eta)], linestyle = 4, color = black, legend = "Le = 2.0"); plots[display]({p5, p6, p7, p8});



(when  b > 6, for the fix vales of rest of the parameters there must exist the dual solution(i.e two solutions). but an error is occure over here. Some one help me please?)

my email address is 
r.haq.qau@gmail.com