Pleas help me to run this simple code

T K
17 Feb 2026
1 Answer
Updated 17 Feb 2026
13 Views (30 days)
Ran in:

0 votes

Ran in:
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

Sign in to answer this question.

Answers (1)

Walter Roberson
Walter Roberson on 17 Feb 2026
Moved: Walter Roberson on 17 Feb 2026
Ran in:

0 votes

Ran in:
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

3 Comments
Show 1 older comment Hide 1 older comment

T K
T K on 17 Feb 2026
Edited: T K on 17 Feb 2026
@ 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
T K on 17 Feb 2026
What do you mean component #7 ?
Walter Roberson
Walter Roberson on 17 Feb 2026
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

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Tags

Asked:

T K
T K
on 17 Feb 2026

Commented:

Walter Roberson
Walter Roberson
on 17 Feb 2026

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!