Unable to solve the collocation equations -- a singular Jacobian encountered.

hello every body, I'm facing

           Error using bvp4c (line 251)
           Unable to solve the collocation equations -- a singular Jacobian encountered.
           Error in num_opt_ctrl (line 11)
           sol = bvp4c(@ode, @bc, solinit);

but i don't know how to solve it, here's my code:

function num_opt_ctrl
clc
clear all
T = 100000;%Thrust
M = 1000;%mass
g0 = 1.62;%moon gravity
R = 1737000;%moon radius
D = R/2;%altitude
solinit = bvpinit(linspace(0,1),[R; 0; 0; 0; 2; 2; 3; 4; 100]);

sol = bvp4c(@ode, @bc, solinit);
y = sol.y;
time = y(9)*sol.x;
ut = atan(y(7,:)/y(8,:));

figure(1);
plot(time,y([1 2],:)','-'); hold on;
plot(time, ut, 'k:');
axis([0 time(1,end) -1.5 3]);
text(1.3,2.5,'x_1(t)');
text(1.3,.9,'x_2(t)');
text(1.3,-.5,'u(t)');
xlabel('time');
ylabel('states');
title('Numerical solution');
hold off;

ODE's of augmented states

function dydt = ode(t,y)
T = 100000;%Thrust
M = 1000;%mass
g0 = 1.62;%moon gravity
R = 1737000;%moon radius
D = R/2;%altitude
dydt = y(9)*[y(3);
             y(4)/y(1);
             y(4)^2/y(1) - g0*R^2/y(1)^2 + (T/M)*(y(7)/norm(y(7),y(8)));
             -y(3)*y(4)/y(1) + (T/M)*(y(8)/norm(y(7),y(8)));
             y(6)*y(4)/y(1)^2 + y(7)*(y(4)^2/y(1)^2 - 2*g0*R^2/y(1)^3) - y(8)*y(3)*y(4)/y(1)^2;
             0;
             -y(5) + y(8)*y(4)/y(1);
             -y(6)/y(1) - 2*y(7)*y(4)/y(1) + y(8)*y(3)/y(1);
             0];

boundary conditions: x1(0)=R;x2(0)=0, x3(0)=0, x4(0)=0, x1(tf)=R+D, x3(tf)=0, x4(tf)=sqrt(g0*R^2/(R+D)); lam2(tf)=0; 1 + lam1(tf)*x3(tf) + lam2(tf)*x4(tf)/x1(tf) + ... lam3(tf)*(x4(tf)^2/x1(tf)-g0*R^2/x1(tf)^2+(T/M)*(lam3(tf)/norm(lam3(tf),lam4(tf))))... lam4(tf)*(-x3(tf)*x4(tf)/x1(tf)+(T/M)*(lam4(tf)/norm(lam3(tf),lam4(tf))))

function res = bc(ya,yb)
T = 100000;%Thrust
M = 1000;%mass
g0 = 1.62;%moon gravity
R = 1737000;%moon radius
D = R;%altitude
res = [ ya(1) - R; 
        ya(2); 
        ya(3); 
        ya(4); 
        yb(1) - (R+D); 
        yb(3); 
        yb(4) - sqrt(g0*R^2/(R+D)); 
        yb(6);
        1 + yb(5)*yb(3) + yb(6)*yb(4)/yb(1) + ...
        yb(7)*(yb(4)^2/yb(1)-g0*R^2/yb(1)^2+(T/M)*(yb(7)/norm(yb(7),yb(8)))) + ...
        yb(8)*(-yb(3)*yb(4)/yb(1)+(T/M)*(yb(8)/norm(yb(7),yb(8))))];