Nt=0.5; Nb=0.5; Le=2; Pr=1; lambda=1.5; s=1; A=3;
sol = bvpinit(linspace(0,6,10), [0 0 0 0 0 0 0]);
sol1 = bvp4c(@bvpexam2, @bcexam2, sol);
opts = bvpset('stats','off','RelTol',1e-10);
sol = bvpinit(linspace(0,5,10), [-1 0 0 0 0 0 0]);
sol2 = bvp4c(@bvpexam2, @bcexam2_dual, sol,opts);
plot(x1,y1(3,:),'-'); hold on
function res = bcexam2(y0, yinf)
res= [y0(1)-s; y0(2)-lambda; y0(4)-1; y0(6)-1; yinf(2); yinf(4);yinf(6)];
function res = bcexam2_dual(y0, yinf)
res= [y0(1)-s; y0(2)-lambda; y0(4)-1; y0(6)-1; yinf(2); yinf(4);yinf(6)];
function ysol = bvpexam2(x,y)
yy1 = -(A*y(1)*y(3)-A*(y(2))^2)-y(2)-(x/2)*y(3);
yy2 = -Pr*(A*y(1)*y(5)+(x/2)*y(5)+Nb*y(5)*y(7)+Nt*(y(5))^2);
yy3 = (-Le*(A*(y(1)*y(7)+(x/2)*y(7)))-(Nt/Nb)*( -Pr*(A*y(1)*y(5)+Nb*y(5)*y(7)+Nt*(y(5))^2)));
ysol = [y(2); y(3); yy1;y(5);yy2;y(7);yy3];