sol2 = bvp4c(@bvpexam2, @bcexam2, sol);

how to get second solution from this code.I ploted first solution for different values of alpha.

function  Ibrardual
clc
clear all
Nt=0.5; Nb=0.5; Le=2; Pr=1; alpha=-0.2138; s=1; A=1;
%% solution in structure form
%First solution
sol = bvpinit(linspace(0,2,10), [1 0 0 0 0 0 0]);
sol1 = bvp4c(@bvpexam2, @bcexam2, sol);
x1 = sol1.x;
y1=sol.y;
plot(x1,y1(3,:),'b');
hold on
result=(A)^-1/2*y1(3,1)


%% Here I define residual of boundary conditions
    function res = bcexam2(y0, yinf)        
        res= [y0(1)-s; y0(2)-alpha; y0(4)-1; y0(6)-1;
            yinf(2); yinf(4);yinf(6)];
    end
%% First order ODEs are define here
    function ysol = bvpexam2(x,y) 
        
        yy1  = -A*(y(1)*y(3)-(y(2))^2)-y(2);
        yy2  = -Pr*(A*y(1)*y(5)+Nb*y(5)*y(7)+Nt*(y(5))^2);
        yy3  = (-Le*A*(y(1)*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];
        
    end
end

Answer 1

[Torsten]

 Ibrardual()

function  Ibrardual
Nt=0.5; Nb=0.5; Le=2; Pr=1; Alpha=[-0.2138;-0.4]; s=1; A=1;
%% solution in structure form
%First solution
hold on
for i = 1:numel(Alpha)
  alpha = Alpha(i);
  sol = bvpinit(linspace(0,2,10), [1 0 0 0 0 0 0]);
  sol1 = bvp4c(@bvpexam2, @bcexam2, sol);
  x1 = sol1.x;
  y1 = sol1.y;
  plot(x1,y1(3,:),'b');
end
hold off
 %% Here I define residual of boundary conditions
    function res = bcexam2(y0, yinf)        
res= [y0(1)-s; y0(2)-alpha; y0(4)-1; y0(6)-1;
    yinf(2); yinf(4);yinf(6)];
    end
%% First order ODEs are define here
    function ysol = bvpexam2(x,y) 
                 
yy1  = -A*(y(1)*y(3)-(y(2))^2)-y(2);
yy2  = -Pr*(A*y(1)*y(5)+Nb*y(5)*y(7)+Nt*(y(5))^2);
yy3  = (-Le*A*(y(1)*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];
        
    end
end