Pleas help me to run this simple code

proj()

Error using bvp4c (line 196)
Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in solution>proj (line 17)
sol= bvp4c(@projfun,@projbc,solinit,options);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

function sol= proj
myLegend1 = {};
k0=386; ce=3.831*10^2; mu=38.6*10^9;alfat=1.78*10^-5; rho=89.54*10^2; lamda=77.6*10^9;taw=0.5;Tnot=2.93*10^2;
c0=sqrt((lamda+2*mu)/(rho)); Betanot=(3*lamda+2*mu)*alfat; a1=mu/(lamda+2*mu);a2=(mu+lamda)/(lamda+2*mu);a3=(Betanot*Tnot)/(lamda+2*mu);omega=(rho*ce)/(k0);
a4=lamda/(lamda+2*mu);a5=(k0*omega*c0^2)/(k0);a6=(rho*ce*c0^2)/(k0);
a7=(Betanot*c0^2)/(k0); a8=a6*taw; a9=a7*taw; a10=rho*ce*taw*omega*c0^4/(k0); a11=Betanot*taw*omega*c0^4/(k0);w=rho*ce/(k0);
rr = [0.2 0.4 0.6];
for i =1:numel(rr)
    a= rr(i);
    s=1;h=1;
    y0 = [1,0,0,1,1,0,0];
    options =bvpset('stats','on','RelTol',1e-4);
    m = linspace(0,5);
    solinit = bvpinit(m,y0);
    sol= bvp4c(@projfun,@projbc,solinit,options);
    figure(1)
    plot(sol.x,(sol.y(5,:)))
    grid on,hold on
    myLegend1{i}=['alfa= ',num2str(rr(i))];
    
end
figure(1)
legend(myLegend1)
hold on
    function dy = projfun(x,y)
        dy = zeros(7,1);
        v = y(1);
        dv = y(2);
        u = y(3);
        du = y(4);
        t = y(5);
        dt = y(6);
        ddt = y(7);
        dy(1) = dv;
        dy(2)=(1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a*(s^2+a1*h^2)*t-a2*s*h*a*(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h));
        dy(3)=du;
        dy(4) =(1/(a*(s^2+a1*h^2-(x+1+1)^2)*t-(s^2+a1*h^2-(x+1+1)^2)))*((a2*s*h-a*a2*s*h*t)*((1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a*(s^2+a1*h^2)*t-a2*s*h*a*(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h)))+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dv);
        dy(5)=dt;
        dy(6)=ddt;
        dy(7)=(1/(a*(a5*s^2*(x+2)+a5*h^2*(x+2)-a8*(x+2)^3)*t+a5*s^2*(x+2)^2+a5*h^2*(x+2)-a8*(x+2)^3))*(-(2*a7*(x+2)-2*a*a11*(x+2))*((1/(a*(s^2+a1*h^2-(x+1+1)^2)*t-(s^2+a1*h^2-(x+1+1)^2)))*((a2*s*h-a*a2*s*h*t)*((1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a*(s^2+a1*h^2)*t-a2*s*h*a*(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h)))+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dv))+(s^2+h^2+2*a5*t-a6*(x+2)^2)*ddt-(a*(s^2+h^2+4*a5-a6*(x+2)^2))*t*ddt-(a*(s^2+h^2+a5+a5))*dt*dt-a*(x+2)*(a5*s^2+a5*h^2+a10*(x+2)^2)*dt*ddt-(2*a7*(x+2)-2*a*a11*(x+2))*((1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a*(s^2+a1*h^2)*t-a2*s*h*a*(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h)))+2*a*a7*(x+2)*t*((1/(a*(s^2+a1*h^2-(x+1+1)^2)*t-(s^2+a1*h^2-(x+1+1)^2)))*((a2*s*h-a*a2*s*h*t)*((1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a*(s^2+a1*h^2)*t-a2*s*h*a*(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h)))+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dv))+(2*a*a7*(x+2))*t*((1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a*(s^2+a1*h^2)*t-a2*s*h*a*(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h))));
    end
    function res = projbc(ya,yb)
        res = [ya(1);
            ya(3);
            ya(5)-1;
            ya(7);
            yb(1);
            yb(3);
            yb(5);
            ];
    end
end

Answer 1

[Walter Roberson]

Your equations are going to +/- infinity in the last component.

Below, max_dy is printed out each time max(abs(dy)) increases

proj()

max dy increased to -1.03787e+09 at component #7
max dy increased to -1.04008e+09 at component #7
max dy increased to -1.0901e+09 at component #7
max dy increased to -1.17508e+09 at component #7
max dy increased to -1.31614e+09 at component #7
max dy increased to -1.56033e+09 at component #7
max dy increased to -2.03538e+09 at component #7
max dy increased to -3.24771e+09 at component #7
max dy increased to -1.14454e+10 at component #7
max dy increased to 1.51966e+10 at component #7
max dy increased to 1.62624e+10 at component #7
max dy increased to -9.60812e+13 at component #7
max dy increased to -2.15787e+14 at component #7
max dy increased to -3.60165e+14 at component #7
max dy increased to -5.29562e+14 at component #7
max dy increased to -7.23281e+14 at component #7
max dy increased to -9.3935e+14 at component #7
max dy increased to -1.17413e+15 at component #7
max dy increased to -1.42209e+15 at component #7
max dy increased to -1.67524e+15 at component #7
max dy increased to -1.92399e+15 at component #7
max dy increased to -2.14796e+15 at component #7
max dy increased to -2.22608e+15 at component #7
max dy increased to -2.48603e+15 at component #7
max dy increased to -2.64374e+15 at component #7
max dy increased to 3.73381e+15 at component #7
max dy increased to 6.47483e+15 at component #7
max dy increased to 9.99645e+15 at component #7
max dy increased to 1.44491e+16 at component #7
max dy increased to 2.00057e+16 at component #7
max dy increased to 2.68613e+16 at component #7
max dy increased to 3.52405e+16 at component #7
max dy increased to 4.53984e+16 at component #7
max dy increased to 5.76213e+16 at component #7
max dy increased to 7.22314e+16 at component #7
max dy increased to 8.95989e+16 at component #7
max dy increased to 1.10133e+17 at component #7
max dy increased to 1.343e+17 at component #7
max dy increased to 1.62614e+17 at component #7
max dy increased to 1.95671e+17 at component #7
max dy increased to 2.34109e+17 at component #7
max dy increased to 2.78674e+17 at component #7
max dy increased to 3.30185e+17 at component #7
max dy increased to 3.89547e+17 at component #7
max dy increased to 4.57799e+17 at component #7
max dy increased to 5.36075e+17 at component #7
max dy increased to 6.25684e+17 at component #7
max dy increased to 7.28032e+17 at component #7
max dy increased to 8.44769e+17 at component #7
max dy increased to 9.77667e+17 at component #7
max dy increased to 1.12873e+18 at component #7
max dy increased to 1.30023e+18 at component #7
max dy increased to 1.49472e+18 at component #7
max dy increased to 1.71502e+18 at component #7
max dy increased to 1.96437e+18 at component #7
max dy increased to 2.24637e+18 at component #7
max dy increased to 2.56499e+18 at component #7
max dy increased to 2.92487e+18 at component #7
max dy increased to 3.3311e+18 at component #7
max dy increased to 3.78963e+18 at component #7
max dy increased to 4.3069e+18 at component #7
max dy increased to 4.89047e+18 at component #7
max dy increased to 5.54868e+18 at component #7
max dy increased to 6.2911e+18 at component #7
max dy increased to 7.12886e+18 at component #7
max dy increased to 8.07401e+18 at component #7
max dy increased to 9.14104e+18 at component #7
max dy increased to 1.0346e+19 at component #7
max dy increased to 1.17075e+19 at component #7
max dy increased to 1.32462e+19 at component #7
max dy increased to 1.49867e+19 at component #7
max dy increased to 1.69571e+19 at component #7
max dy increased to 1.91888e+19 at component #7
max dy increased to 2.17194e+19 at component #7
max dy increased to 2.45911e+19 at component #7
max dy increased to 2.78533e+19 at component #7
max dy increased to 3.15632e+19 at component #7
max dy increased to 3.57857e+19 at component #7
max dy increased to 4.06003e+19 at component #7
max dy increased to 4.60941e+19 at component #7
max dy increased to 5.2373e+19 at component #7
max dy increased to 5.9558e+19 at component #7
max dy increased to 6.77929e+19 at component #7
max dy increased to 7.72436e+19 at component #7
max dy increased to 8.8113e+19 at component #7
max dy increased to 1.00628e+20 at component #7
max dy increased to 1.15065e+20 at component #7
max dy increased to 1.31758e+20 at component #7
max dy increased to 1.51085e+20 at component #7
max dy increased to 1.73516e+20 at component #7
max dy increased to 1.99599e+20 at component #7
max dy increased to 2.30006e+20 at component #7
max dy increased to 2.65529e+20 at component #7
max dy increased to 3.07146e+20 at component #7
max dy increased to 3.56019e+20 at component #7
max dy increased to 4.13585e+20 at component #7
max dy increased to 4.81598e+20 at component #7
max dy increased to 5.62216e+20 at component #7
max dy increased to 6.58104e+20 at component #7
max dy increased to 7.72611e+20 at component #7
max dy increased to 9.09916e+20 at component #7
max dy increased to 1.07528e+21 at component #7
max dy increased to 1.27549e+21 at component #7
max dy increased to 1.51921e+21 at component #7
max dy increased to 4.7741e+22 at component #7
max dy increased to 6.71942e+22 at component #7
max dy increased to 9.10873e+22 at component #7
max dy increased to 1.19593e+23 at component #7
max dy increased to 1.52414e+23 at component #7
max dy increased to 1.88675e+23 at component #7
max dy increased to 2.2652e+23 at component #7
max dy increased to 2.62798e+23 at component #7
max dy increased to 2.91515e+23 at component #7
max dy increased to 2.99859e+23 at component #7
max dy increased to -5.75886e+23 at component #7
max dy increased to -1.14348e+24 at component #7
max dy increased to -1.95779e+24 at component #7
max dy increased to -3.09227e+24 at component #7
max dy increased to -4.63926e+24 at component #7
max dy increased to -6.71261e+24 at component #7
max dy increased to -9.45247e+24 at component #7
max dy increased to -1.30304e+25 at component #7
max dy increased to -1.76537e+25 at component #7
max dy increased to -2.35782e+25 at component #7
max dy increased to -3.10957e+25 at component #7
max dy increased to -4.05812e+25 at component #7
max dy increased to -5.24738e+25 at component #7
max dy increased to -6.72817e+25 at component #7
max dy increased to -8.56113e+25 at component #7
max dy increased to -1.08236e+26 at component #7
max dy increased to -1.35987e+26 at component #7
max dy increased to -1.69926e+26 at component #7
max dy increased to -2.1122e+26 at component #7
max dy increased to -2.61423e+26 at component #7
max dy increased to -3.22037e+26 at component #7
max dy increased to -3.95277e+26 at component #7
max dy increased to -4.83392e+26 at component #7
max dy increased to -5.89055e+26 at component #7
max dy increased to -7.15744e+26 at component #7
max dy increased to -8.67043e+26 at component #7
max dy increased to -1.04794e+27 at component #7
max dy increased to -1.2631e+27 at component #7
max dy increased to -1.51972e+27 at component #7
max dy increased to -1.82441e+27 at component #7
max dy increased to -2.18602e+27 at component #7
max dy increased to -2.61536e+27 at component #7
max dy increased to -3.12451e+27 at component #7
max dy increased to -3.72798e+27 at component #7
max dy increased to -4.44332e+27 at component #7
max dy increased to -5.2913e+27 at component #7
max dy increased to -6.29377e+27 at component #7
max dy increased to -7.48353e+27 at component #7
max dy increased to -8.89177e+27 at component #7
max dy increased to -1.05643e+28 at component #7
max dy increased to -1.2545e+28 at component #7
max dy increased to -1.48981e+28 at component #7
max dy increased to -1.76901e+28 at component #7
max dy increased to -2.10081e+28 at component #7
max dy increased to -2.4962e+28 at component #7
max dy increased to -2.96602e+28 at component #7
max dy increased to -3.52729e+28 at component #7
max dy increased to -4.19687e+28 at component #7
max dy increased to -4.99824e+28 at component #7
max dy increased to -5.95495e+28 at component #7
max dy increased to -7.10367e+28 at component #7
max dy increased to -8.48394e+28 at component #7
max dy increased to -1.01411e+29 at component #7
max dy increased to -1.21408e+29 at component #7
max dy increased to -1.45518e+29 at component #7
max dy increased to -1.74677e+29 at component #7
max dy increased to -2.10016e+29 at component #7
max dy increased to -2.52833e+29 at component #7
max dy increased to -3.05063e+29 at component #7
max dy increased to -3.68578e+29 at component #7
max dy increased to -4.46245e+29 at component #7
max dy increased to -5.4126e+29 at component #7
max dy increased to -6.57881e+29 at component #7
max dy increased to -8.01162e+29 at component #7
max dy increased to -9.78251e+29 at component #7
max dy increased to -1.19656e+30 at component #7
max dy increased to -1.46725e+30 at component #7
max dy increased to -1.80408e+30 at component #7
max dy increased to -2.22272e+30 at component #7
max dy increased to -2.74642e+30 at component #7
max dy increased to -3.40196e+30 at component #7
max dy increased to -4.22666e+30 at component #7
max dy increased to -5.26523e+30 at component #7
max dy increased to -6.58025e+30 at component #7
max dy increased to -8.24736e+30 at component #7
max dy increased to -1.03717e+31 at component #7
max dy increased to -1.3088e+31 at component #7
max dy increased to -1.65747e+31 at component #7
max dy increased to -2.10686e+31 at component #7
max dy increased to -2.68927e+31 at component #7
max dy increased to -3.44745e+31 at component #7
max dy increased to -4.43946e+31 at component #7
max dy increased to -5.74656e+31 at component #7
max dy increased to -7.47928e+31 at component #7
max dy increased to -7.14289e+46 at component #7
max dy increased to -4.43626e+54 at component #7
max dy increased to -5.19348e+54 at component #7
max dy increased to 3.19614e+63 at component #7
max dy increased to 3.56268e+75 at component #7
max dy increased to -1.02917e+79 at component #7
max dy increased to -9.74423e+83 at component #7
max dy increased to -6.39195e+86 at component #7
max dy increased to -1.08617e+96 at component #7
max dy increased to -2.50327e+139 at component #7
max dy increased to -2.5451e+147 at component #7
max dy increased to -1.82551e+148 at component #7
max dy increased to -1.82551e+148 at component #7
max dy increased to -4.68542e+152 at component #7
max dy increased to 6.2098e+163 at component #7
max dy increased to 2.86595e+179 at component #7
max dy increased to 7.25633e+239 at component #7
max dy increased to -Inf at component #2

Error using bvp4c (line 196)
Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in solution>proj (line 17)
sol= bvp4c(@projfun,@projbc,solinit,options);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

function sol= proj
myLegend1 = {};
k0=386; ce=3.831*10^2; mu=38.6*10^9;alfat=1.78*10^-5; rho=89.54*10^2; lamda=77.6*10^9;taw=0.5;Tnot=2.93*10^2;
c0=sqrt((lamda+2*mu)/(rho)); Betanot=(3*lamda+2*mu)*alfat; a1=mu/(lamda+2*mu);a2=(mu+lamda)/(lamda+2*mu);a3=(Betanot*Tnot)/(lamda+2*mu);omega=(rho*ce)/(k0);
a4=lamda/(lamda+2*mu);a5=(k0*omega*c0^2)/(k0);a6=(rho*ce*c0^2)/(k0);
a7=(Betanot*c0^2)/(k0); a8=a6*taw; a9=a7*taw; a10=rho*ce*taw*omega*c0^4/(k0); a11=Betanot*taw*omega*c0^4/(k0);w=rho*ce/(k0);
rr = [0.2 0.4 0.6];
for i =1:numel(rr)
    a= rr(i);
    s=1;h=1;
    y0 = [1,0,0,1,1,0,0];
    options =bvpset('stats','on','RelTol',1e-4);
    m = linspace(0,5);
    solinit = bvpinit(m,y0);
    sol= bvp4c(@projfun,@projbc,solinit,options);
    figure(1)
    plot(sol.x,(sol.y(5,:)))
    grid on,hold on
    myLegend1{i}=['alfa= ',num2str(rr(i))];
    
end
figure(1)
legend(myLegend1)
hold on
    function dy = projfun(x,y)
        global max_dy
        if isempty(max_dy); max_dy = 0; end
        dy = zeros(7,1);
        v = y(1);
        dv = y(2);
        u = y(3);
        du = y(4);
        t = y(5);
        dt = y(6);
        ddt = y(7);
        dy(1) = dv;
        dy(2)=(1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a*(s^2+a1*h^2)*t-a2*s*h*a*(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h));
        dy(3)=du;
        dy(4) =(1/(a*(s^2+a1*h^2-(x+1+1)^2)*t-(s^2+a1*h^2-(x+1+1)^2)))*((a2*s*h-a*a2*s*h*t)*((1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a*(s^2+a1*h^2)*t-a2*s*h*a*(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h)))+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dv);
        dy(5)=dt;
        dy(6)=ddt;
        dy(7)=(1/(a*(a5*s^2*(x+2)+a5*h^2*(x+2)-a8*(x+2)^3)*t+a5*s^2*(x+2)^2+a5*h^2*(x+2)-a8*(x+2)^3))*(-(2*a7*(x+2)-2*a*a11*(x+2))*((1/(a*(s^2+a1*h^2-(x+1+1)^2)*t-(s^2+a1*h^2-(x+1+1)^2)))*((a2*s*h-a*a2*s*h*t)*((1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a*(s^2+a1*h^2)*t-a2*s*h*a*(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h)))+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dv))+(s^2+h^2+2*a5*t-a6*(x+2)^2)*ddt-(a*(s^2+h^2+4*a5-a6*(x+2)^2))*t*ddt-(a*(s^2+h^2+a5+a5))*dt*dt-a*(x+2)*(a5*s^2+a5*h^2+a10*(x+2)^2)*dt*ddt-(2*a7*(x+2)-2*a*a11*(x+2))*((1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a*(s^2+a1*h^2)*t-a2*s*h*a*(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h)))+2*a*a7*(x+2)*t*((1/(a*(s^2+a1*h^2-(x+1+1)^2)*t-(s^2+a1*h^2-(x+1+1)^2)))*((a2*s*h-a*a2*s*h*t)*((1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a*(s^2+a1*h^2)*t-a2*s*h*a*(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h)))+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dv))+(2*a*a7*(x+2))*t*((1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a*(s^2+a1*h^2)*t-a2*s*h*a*(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h))));
        [~, pos] = max(abs(dy));
        T = dy(pos(1));
        if abs(T) > abs(max_dy)
            max_dy = T;
            fprintf('max dy increased to %g at component #%d\n', max_dy, pos(1));
        end
    end
    function res = projbc(ya,yb)
        res = [ya(1);
            ya(3);
            ya(5)-1;
            ya(7);
            yb(1);
            yb(3);
            yb(5);
            ];
    end
end

[T K]

@ Walter Roberson. Thank you for help me. This error is eliminated by the disappearance of terms for the component that causes infinity in the description of the latter equation?

[T K]

What do you mean component #7 ?

[Walter Roberson]

Component #7 refers to dy(7)

Component #2 refers to dy(2)

"This error is eliminated by the disappearance of terms for the component that causes infinity in the description of the latter equation?"

No, not likely. If you remove those components, then you will be working with a model that does not match the equations.

You need to verify that the implementing equations are correct.

You could consider initially implementing the equations using the Symbolic Toolbox, and then converting the equations to numeric form using odeFunction