IVP not Lipschitz --> BVP

In some problems a first-order ODE IVP does not pose a unique solution because it is not Lipschitz continuous. Instead it is proposed that one introduces a second boundary condition such the problem is a BVP.

Could anyone help in guiding me to the structure of such a code in Matlab.

If I take for the simplicity (just ignore the specifics, more interested in how it is setting up) dy/dx = y^2 + y*x*a with boundary condition Y(0)=c Y(5)=d

Would is wrong with writing it up like this - is there a more simple way of having first-order ODE with two boundaries?

xlow = 0; xhigh = 5; solinit = bvpinit(linspace(xlow,xhigh,10),[0,1]); sol = bvp4c(@(x,y)[y(1)^2 + y*x*a],@(ya,yb)[ya(1)-c; yb(1)-d],solinit); x = linspace(xlow, xhigh, 100); y = deval(sol, x); plot(x, y(1,:))

Thanks.