@acer 

equ:

Graph:

code:

restart; with(plots); fcns := {C(eta), T(eta), f(eta)}; bet := 1.0; M := 3.0; pr := 6.0; nb := .2; nt := .7; le := 1.0; bi := 3.0; S := 0.; lam := 1.0; N := 10; sys := (1+1/bet)*(diff(f(eta), eta, eta, eta))-(M*M)*(diff(f(eta), eta))-2*(diff(f(eta), eta))*(diff(f(eta), eta))+f(eta)*(diff(f(eta), eta, eta)) = 0, diff(T(eta), eta, eta)+pr*(f(eta)*(diff(T(eta), eta))-(diff(T(eta), eta))*T(eta))+pr*nb*(diff(T(eta), eta))*(diff(C(eta), eta))+pr*nt*(diff(T(eta), eta))*(diff(T(eta), eta)) = 0, diff(C(eta), eta, eta)+le*pr*(f(eta)*(diff(C(eta), eta))-(diff(f(eta), eta))*C(eta))+nt*(diff(T(eta), eta, eta))/nb = 0; bc := f(0) = S, (D(f))(0) = lam, (D(f))(N) = 0, (D(T))(0) = -bi*(1-T(0)), T(N) = 0, C(0) = 1, C(N) = 0; R := dsolve(eval({bc, sys}), fcns, type = numeric, method = bvp[midrich], output = operator); psi := unapply(eval(x*f(eta)+T(eta), R), x, eta); contourplot(psi(eta, x), eta = 0 .. N, x = 0 .. 5, filledregions = true)

I need the above graph as it is. Have a good day.

@KitonumI receive like this

@Christopher2222  yes exactly what you say, velocity and temperature fields. 

@acer Sorry, I have no idea what u are asking, anyway Thanks for your cooperation. have a good day.

@acer Here U and W both are velocities. y is the eta (eg problem). Here want to plot U wrt y