I am using bvp4c from Matlab to solve a boundary values ODEs problem.
Given the ODEs and boundary conditions, is there any way to have more information on the solutions? How many do I expect? 0, 1, several?
What am I supposed to understand from the bvp4c error message "singular Jacobian encountered"? No solutions at all? Initial guess doesn't fit?
A singular Jacobian means that there is a linear system in there somewhere (probably within a Newton solve) that has multiple (or zero) solutions. This issue may correspond to different solutions having different amounts of material; you might be able to resolve it by adding an equation that governs the amount of material present.
EDIT: MATLAB's bvp4c uses a three-stage Lobatto IIIa formula that is fourth-order accurate. Essentially, there's some default discretization that defines the collocation points of the method, and then the Lobatto formula defines the fourth-order interpolating polynomial on that discretization. Since Lobatto IIIa methods are implicit methods, they necessarily involve solving a nonlinear equation, and this nonlinear equation solve is usually done through a variant of Newton's method, which involves solving a linear system at each iteration.
Without knowing the details of the BVP solver algorithm at the time of my original response, the existence of the linear solve in the algorithm could be inferred because:
EDIT 2: There's no general theory for existence and uniqueness of two-point boundary value problems (TPBVP). There are results for specific cases.
In terms of intuition, if you have a TPBVP expressed in first-order form with $n$ equations, you need $n$ linearly independent boundary conditions. A common case involving multiple solutions occurs when you have all Neumann boundary conditions (i.e., all boundary conditions involve derivatives only). Usually, if you have at least one Dirichlet condition, it suggests there might be a solution to your problem (although it may not be unique).
As I wrote above, if you get an error message "singular Jacobian encountered", one common situation is that there should be an additional constraint on your system in the form of another equation for your boundary value problem; adding additional constraint equations is a common way of determining a unique solution to a problem with only derivative boundary conditions, and these constraints could be things like conservation of mass (so the integral of the solution over the domain must be a constant).
It could be a problem with your initial guess, too, but without knowing more about intermediate outputs from your problem or the specific problem you are solving, it is difficult to diagnose the precise cause. Generally, if there is a problem with your initial guess, when you look at the iterates for the nonlinear solves in the algorithm, they fail to converge and the residual function values stagnate; it's not clear to me from looking at the documentation if there are MATLAB options to output this information. You could also have a singular TPBVP, and you may have to resort to special methods to solve your problem. I find other software gives you more control over the numerical methods and access to solver information, but it generally comes at the cost of having to implement more of the numerical methods yourself. You could use PETSc, for instance, which would give you access to initial-value problem ODE solvers, nonlinear solvers, and linear solvers, along with tons of solver information that would help you better diagnose these problems, but then you would have to write your own boundary value problem solver (I'd suggest a multiple shooting method).
