function main
format('long');
gg=['r','k','b','g','m','c','y','r.','m.','k.'];
Pr=10; Le=10;
Nb=.1;
xa=0;xb=2;
Nt=1e-3;
Nt=[0.1:0.1:0.5];
solinit=bvpinit(linspace(xa,xb,1000),[0 1 1 1 0 0 0]);
sol=bvp4c(@ode,@bc,solinit);
xint=linspace(xa,xb,100);
sxint=deval(sol,xint);
function res=bc(ya,yb)
res=[ya(1); ya(2)-1; ya(4)-1; ya(6)-1; yb(2); yb(4); yb(6)];
end
function dydx=ode(x,y)
dydx=[y(2); y(3); y(2)^2-y(3)*y(1); y(5); -Pr*(y(1)*y(5)+Nb*y(5)*y(7)+Nt*y(5)^2);
y(7); -Le*y(7)*y(1)+(Nt/Nb)*Pr*(y(1)*y(5)+Nb*y(5)*y(7)+Nt*y(5)^2)];
end
T0 = deval(sol,0);
-T0(5);
plot(Nt,-T0(5))