function answers_5_14_2019
dydx=@(x,y) f(x,y);
Uh=3;
bcFunc=@(ya,yb) bc(ya,yb,Uh);
yinit=@(x)[x; 1];
L=10;
solinit=bvpinit(linspace(0,L,20),yinit);
sol=bvp4c(dydx,bcFunc,solinit);
figure; plot(sol.x,sol.y(1,:)); grid;
end

function dydx = f(x,y)
a = .1; b = 1; c = 0; d = 1e-3*x; e = 0; f = 0;
  dydx = [y(2);
          -1/(b*y(2))*(a*y(2)^2 + c*y(1) +d*y(2) + e*y(1)^2 + f)];
end
% -----------------------------------------------------------------------
function res = bc(ya, yb, Uh)
% Boundary conditions
  res = [ya(1); yb(1)-Uh];
end