hold on
for k = 1:4;
Ha = k;
init = bvpinit(linspace(-1,1,20),[0 0]);
sol=bvp4c(@(x,y)rhs_bvp(x,y,Ha),@bc_bvp,init);
x = linspace(-1,1,20);
BS(:,:,k) =deval(sol,x);
plot(BS(1,:,k),x, '--');
xlim([-1 1])
end
hold off
function rhs=rhs_bvp(x,y,Ha)
rhs = [y(2); ((Ha^2+(1/3))/(1+0.3))*y(1)];
end
function bc=bc_bvp(yl,yr)
bc = [yl(1)-(0.1*1.3)*yl(2)-1; yr(1)+(0.1*1.3)*yr(2)+1];
end
