Question: how to plot f''(0), theta'(0) and phi'(0) against two parameters

restart; b := 6; de1 := diff(f(eta), eta, eta, eta)+f(eta)*(diff(f(eta), eta, eta))-M^2*(diff(f(eta), eta))-(diff(f(eta), eta))^2 = 0, (1+(4/3)*Nr)*(diff(theta(eta), eta, eta))+epsilon*(theta(eta)*(diff(theta(eta), eta, eta))+epsilon*(diff(theta(eta), eta))^2)+Pr*(s*theta(eta)+f(eta)*(diff(theta(eta), eta))+Nb*(diff(theta(eta), eta))*(diff(phi(eta), eta))+Nt*(diff(theta(eta), eta))^2) = 0, diff(phi(eta), eta, eta)+Le*f(eta)*(diff(phi(eta), eta))+Nb*(diff(theta(eta), eta, eta))/Nt = 0, f(0) = 0, (D(f))(0) = 1, (D(theta))(0) = -gamma*(1+theta(0)), phi(0) = 1, (D(f))(b) = 0, theta(b) = 0, phi(b) = 0; d1 := subs(epsilon = 0.5e-1, s = 0.5e-1, gamma = 0.5e-1, Nr = .3, Pr = .7, Nb = .2, Nt = .1, M = .5, Le = 1, [de1]); d2 := subs(epsilon = 0.5e-1, s = 0.5e-1, gamma = 0.5e-1, Nr = .3, Pr = .7, Nb = .2, Nt = .1, M = 1, Le = 1, [de1]); d3 := subs(epsilon = 0.5e-1, s = 0.5e-1, gamma = 0.5e-1, Nr = .3, Pr = .7, Nb = .2, Nt = .1, M = 1.5, Le = 1, [de1]); d4 := subs(epsilon = 0.5e-1, s = 0.5e-1, gamma = 0.5e-1, Nr = .3, Pr = .7, Nb = .2, Nt = .1, M = 2, Le = 1, [de1]); da1 := dsolve(d1, numeric, output = listprocedure); da2 := dsolve(d2, numeric, output = listprocedure); da3 := dsolve(d3, numeric, output = listprocedure); da4 := dsolve(d4, numeric, output = listprocedure); da1(0); da2(0); da3(0); da4(0); with(plots); p1 := odeplot(da1, [eta, diff(f(eta), eta)], linestyle = 1, color = black, legend = "M=0.5"); p2 := odeplot(da2, [eta, diff(f(eta), eta)], linestyle = 2, color = black, legend = "M=1.0"); p3 := odeplot(da3, [eta, diff(f(eta), eta)], linestyle = 3, color = black, legend = "M=1.5"); p4 := odeplot(da4, [eta, diff(f(eta), eta)], linestyle = 4, color = black, legend = "M=2.0"); plots[display]({p1, p2, p3, p4}); p5 := odeplot(da1, [eta, theta(eta)], linestyle = 1, color = black, legend = "M=0.5"); p6 := odeplot(da2, [eta, theta(eta)], linestyle = 2, color = black, legend = "M=1.0"); p7 := odeplot(da3, [eta, theta(eta)], linestyle = 3, color = black, legend = "M=1.5"); p8 := odeplot(da4, [eta, theta(eta)], linestyle = 4, color = black, legend = "M=2.0"); plots[display]({p5, p6, p7, p8}); p9 := odeplot(da1, [eta, phi(eta)], linestyle = 1, color = black, legend = "M=0.5"); p10 := odeplot(da2, [eta, phi(eta)], linestyle = 2, color = black, legend = "M=1.0"); p11 := odeplot(da3, [eta, phi(eta)], linestyle = 3, color = black, legend = "M=1.5"); p12 := odeplot(da4, [eta, phi(eta)], linestyle = 4, color = black, legend = "M=2.0"); plots[display]({p10, p11, p12, p9});

Kindly find this article from Google. It is easily available see its fig 6, 7 and 8. 
here,Nux/Rax1/4 =theta ‘(0)
Shx/Rax1/4 =phi ‘(0)

Natural convection flow of a nanofluid over a vertical plate with uniform surface heat flux

 how I will draw the graph for f’’(0) against Pr   for Nr = 0.1, 0.3, 0.5  taking rest of the parameter fix

similarly for 
for  theta’(0) against Pr   for Nr = 0.1, 0.3, 0.5  taking rest of the parameter fix
for  phi’(0)    against Pr   for Nr = 0.1, 0.3, 0.5  taking rest of the parameter fix
  for above mention code.