Pleae i can solve this error broblem

Error in proj/projfun (line 48)

dy(2)=(1/((((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t)*(((s^2+a1*h^2-(x+b+1)*(x+b+1)))-a*(s^2+a1*h^2-(x+b+1)*(x+b+1))*t))-(a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h(a4+a1)*dt*dE)-(((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*s*h*(a4+a1)*dt*du));

proj()

Array indices must be positive integers or logical values.

Error in solution>proj/projfun (line 50)
dy(2)=(1/((((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t)*(((s^2+a1*h^2-(x+b+1)*(x+b+1)))-a*(s^2+a1*h^2-(x+b+1)*(x+b+1))*t))-(a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h(a4+a1)*dt*dE)-(((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*s*h*(a4+a1)*dt*du));
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in bvparguments (line 96)
testODE = ode(x1,y1,odeExtras{:});
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in bvp4c (line 119)
bvparguments(solver_name,ode,bc,solinit,options,varargin);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in solution>proj (line 20)
sol= bvp4c(@projfun,@projbc,solinit,options);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

function sol= proj
myLegend1 = {};
myLegend2 = {};
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.4 0.5];
for i =1:numel(rr)
    a= rr(i);
    s=0.0001;h=0.0001;b=0.1;
    y0 = [0,0,0,0,0,0,0];
    options =bvpset('stats','on','RelTol',1e-4);
    m = linspace(0,1);
    
    solinit = bvpinit(m,y0);
    sol= bvp4c(@projfun,@projbc,solinit,options);
    
    figure(1)
    plot(sol.x,abs(sol.y(5,:)))
    grid on,hold on
    myLegend1{i}=['alfa= ',num2str(rr(i))];
    figure(2)
    plot(sol.x,(sol.y(5,:)))
    title('Temperature')
    grid on,hold on
    myLegend2{i}=['alfa= ',num2str(rr(i))];
    i=i+1;
    
end
figure(1)
legend(myLegend1)
hold on
figure(2)
legend(myLegend2)

    function dy = projfun(x,y)
        dy = zeros(7,1);
        E  = y(1);
        dE = y(2);
        u  = y(3);
        du = y(4);
        t  = y(5);
        dt = y(6);
        ddt=y(7);
        dy(1) = dE;
        dy(2)=(1/((((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t)*(((s^2+a1*h^2-(x+b+1)*(x+b+1)))-a*(s^2+a1*h^2-(x+b+1)*(x+b+1))*t))-(a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h(a4+a1)*dt*dE)-(((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*s*h*(a4+a1)*dt*du));
        dy(3)=du;
        dy(4) =(1/((s^2+a1*h^2-(x+b+1).^2)-(a*(s^2+a1*h^2-(x+b+1).^2))*t))*((a2*s*h-(a*a2*s*h)*t)*((1/((((s^2+a1*h^2-(x+b+1).^2)-(a*(s^2+a1*h^2-(x+b+1).^2))*t)*(((s^2+a1*h^2-(x+b+1).^2))-a*(s^2+a1*h^2-(x+b+1).^2)*t))-(a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h(a4+a1)*dt*dE)-(((s^2+a1*h^2-(x+b+1).^2)-(a*(s^2+a1*h^2-(x+b+1).^2))*t))*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*s*h*(a4+a1)*dt*du)))+(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE));
        
        dy(5)=dt;
        dy(6)=ddt;
        dy(7)=(1/((a*(x+b+1)*(a5*(s^2+h^2)-a8*(x+b+1).^2))*t-((x+b+1).*(a5*s^2+a5*h^2-a8*(x+b+1).^2))))*((s^2+h^2+2*a5*s*+2*a5*h-a6*(x+b+1).^2)*ddt+((2*a*a7*(x+b+1))*t-2*(x+b+1)*(a7-a*a11))*((1/((((s^2+a1*h^2-(x+b+1).^2)-(a*(s^2+a1*h^2-(x+b+1).^2))*t)*(((s^2+a1*h^2-(x+b+1).^2))-a*(s^2+a1*h^2-(x+b+1).^2)*t))-(a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h(a4+a1)*dt*dE)-(((s^2+a1*h^2-(x+b+1).^2)-(a*(s^2+a1*h^2-(x+b+1).^2))*t))*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*s*h*(a4+a1)*dt*du)))+((a*2*a7)*t-(2*(x+b+1).*(a7-2*a*a11)))*((1/((s^2+a1*h^2-(x+b+1).^2)-(a*(s^2+a1*h^2-(x+b+1).^2))*t))*((a2*s*h-(a*a2*s*h)*t)*((1/((((s^2+a1*h^2-(x+b+1).^2)-(a*(s^2+a1*h^2-(x+b+1).^2))*t)*(((s^2+a1*h^2-(x+b+1).^2))-a*(s^2+a1*h^2-(x+b+1).^2)*t))-(a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h(a4+a1)*dt*dE)-(((s^2+a1*h^2-(x+b+1).^2)-(a*(s^2+a1*h^2-(x+b+1).^2))*t))*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*s*h*(a4+a1)*dt*du)))+(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)))-a*(s^2+h^2+2*a5*s+2*a5*h-a6*(x+b+1).^2)*t*ddt-a*(s^2+h^2+a5*s+a5*h)*dt*dt-a*(x+b+1)*(a5*s^2+a5*h^2+a10*(x+b+1).^2)*dt*ddt);
        
        
    end
    function res = projbc(ya,yb)
        res = [ya(1)-1;
            ya(3)-1;
            ya(5)-1;
            ya(7);
            yb(1);
            yb(3);
            yb(5);];
    end
end

Answer 1

[Walter Roberson]

Your code for dy(2) includes the sub-expression

h(a4+a1)

which is a request to index h at the location designated by a4+a1

Reminder: MATLAB has absolutely no implied multiplication anywhere (except for the formation of complex constants indicated by a literal number immediately followed by i )

[Torsten]

Same expression h(a4+a1) in dy(2) (once) and dy(7) (twice).

After correction, bvp4c will error with the already known message about the collocation equations.

Answer 2

[Jan]

Inside the code line 48 you find this part:

a*s*h(a4+a1)*dt*dE
%    ^ here

There is a missing operator after the variable h. Then the expression in parenthesis are interpreted as index.

If you split this huge line into readable parts with less than 80 characters, Matlab shows this failing line more precisely:

dy(2)=(1/((((s^2+a1*h^2-(x+b+1)* ...
           (x+b+1))-(a*(s^2+a1*h^2-(x+b+1)* ...
           (x+b+1)))*t)*(((s^2+a1*h^2-(x+b+1)*(x+b+1))) - ...
           a*(s^2+a1*h^2-(x+b+1)*(x+b+1))*t)) - ...
           (a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t) * ...
           (-a3*s*dt+2*a*a3*s*t*dt - a*(s^2+a1*h^2)*dt*du - ...
           a*s*h(a4+a1)*dt*dE) - ...    % ERROR here
           (((s^2+a1*h^2-(x+b+1)*(x+b+1)) - ...
           (a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t)) * ...
           (-a3*h*dt+2*a*a3*h*t*dt - ...
           a*(h^2+a1*s^2) * dt*dE-a*s*h*(a4+a1)*dt*du));

By the way, mu=38.6*10^9 is an expensive power operation, while mu=38.6e9 is a cheap constant.