This is an old question (2014), but it represents an issue that probably continues to puzzle others (including me today).
The documentation is missing an important clue, in that bvp functions need to be able to remeshing the independent variable x. Thus, one cannot use indexing to relate x to the deriative of the dependent variable, (y, dy, d2y, etc.). This rule holds for the initial guess, and for the ode function (containing the system of differential equations) and the residual function (containing the boundary conditions) as well.
That said, one can use anonymous functions to allow for general relationships between x and the y derivatives (see the bvpinit documentation for details). Let me illustrate with an example:
Instead of:
init = bvpint(X,V),
try:
%row vector xV0, and matrix V0 with column vectors for each element of xV0
init = bvpint(X, @(x) interp1(xV0, V0, x))
@(x) indicates that this is anonymous function. interp1 will provide interpolates for the column vector of derivatives at the specified x.
Similar arrangements are needed for handling calculations involving x in the ode and and residual functions.
