The boundary value problem is$$\begin{cases}\dot{x}_i = \begin{cases}(0.5D^{-1}\psi)_i, \text{ if }(0.5D^{-1}\psi)_i \le 0 \\ 0 \text{, otherwise}\end{cases}\\\dot{\psi} = 2\Sigma x \\x(0) = x^0 \\x(T) = 0\end{cases}$$Here $x^0>0$, $T$ is given, $D>0$ is a diagonal matrix, $\Sigma>0$ is a symmetric matrix.  Actually the constraint is $x(t)\ge 0\forall t$, equivalent to $\dot{x}(t) \le 0$.For a 2-dimension case, I try to solve it in MATLAB with this codefunction [tint, traj] = liq2mon_v3(x0, D, sigma, T)  tlow = 0;  thigh = T;  L = length(x0); % L=2 in this example  nsegm = 100; % I know it is for final calc only, not for the mesh  solinit = bvpinit(linspace(tlow, thigh, nsegm), @guess);  sol = bvp5c(@bvp4ode, @bvp4bc, solinit);  tint = linspace(tlow, thigh);  z = deval(sol, tint);  traj = z(1:L, :);  function xpsi = guess(t)    xpsi = zeros(1, L*2);    xpsi(1:L) = (x0*(1.0 - 1.0/T*t))';  end  function dydt = bvp4ode(t, y)    xcomp = y(1:L); % L=2 in this example    psicomp = y(L+1:end);    u    = D \ psicomp / 2.0;    dx   = min(u, 0);     dpsi = 2.0 * sigma * xcomp;    dydt = [dx' dpsi'];   end  function res = bvp4bc(ya, yb)    xcompa = ya(1:2);    xcompb = yb(1:2);    res = [(xcompa - x0)' (xcompb)'];  endendand I get an error *Unable to solve the collocation equations -- a singular Jacobian  encounteredError in ==> liq2mon_v3 at 13sol = bvp5c(@bvp4ode, @bvp4bc, solinit)*However if I change the line dx = min(u, 0);to u(u>0) = 1e-13;dx = u;I do get a solution, with $\min(\min(x)) \approx -8\times 10^{-14}$ for the same $x^0$.Questions:Is there a way to solve this system in MATLAB (or another app) so that the constraint is not violated?Maybe this could be solved in Python, even with the current constraint violation?Why is the Jacobian singular somewhere within MATLAB solver in this case?
Solving constrained BVP, singular Jacobian
matlab;python;ode;boundary conditions;constrained optimization