Please help me to run this simple code

proj()

rr = 1×3

0 0.3000 0.6000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Array indices must be positive integers or logical values.

Error in solution>proj/projfun (line 48)
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));
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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 19)
sol= bvp4c(@projfun,@projbc,solinit,options);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

function sol= proj
clc;clf;clear;
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 0.3 0.6 ]
for i =1:numel(rr)
    a= rr(i); 
    s=0.5;h=0.1;
    y0 = [1,0,0,0,1,0,0];
    options =bvpset('stats','on','RelTol',1e-4);
    m = linspace(0,10);
    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))];
    figure(2)
    plot(sol.x,(sol.y(6,:)))
    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)
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-(2*x+2+1).^2))-(s^2+a1*h^2-(2*x+2+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)))+(((a*2*(x+1)*t)-2*(x+1)-(a*s^2+a*a1*h^2)*t)*du)+(2*a*a3*s*t-a3*s)*dt-(a*a4*s*h+a*a1*s*h)*dt*dv);
        dy(5)=dt;
        dy(6)=ddt;
        dy(7)=(1/(a*(a5*s^2*2*(x+1+1)+a5*h^2*2*(x+1+1))*t-(a5*s^2*2*(x+1+1)+a5*h^2*2*(x+1)-a8*(2*(x+1+1))^3)))*((s^2+h^2+2*a5*h*2-a6*(2*(x+1)).^2-3*a8*2*(x+1+1)*2*(x+1))*ddt);
    end
end
function res= projbc(ya,yb)
res= [ya(1);
    ya(3);
    ya(5)-1;
    ya(7);
    yb(1);
    yb(3);
    yb(5);
    ];
end

[Torsten]

My guess is that there are one or more places where you forgot a multiplication sign, i.e. where you wrote a(...) instead of a*(...). Search for these places.

[dpb]

@T K

   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));

Format your code by wrapping long lines and appropriate spacing, etc., so it's possible to read...there's almost no chance of being able to parse those lines as written without missing something...making shorter, temporary variables wouldn't necessarily be a bad idea, either as one way to reduce the apparent complexity.

[Stephen23]

And get rid of cargo-cult code like this: clc;clf;clear;

What do you think clearing an empty function workspace will do ?

Answer 1

[Stephen23]

You're missing the multiplication operator * in several places. Here's the corrected code with the missing * operators added:

proj()

The solution was obtained on a mesh of 100 points.
The maximum residual is  4.453e-06. 
There were 1303 calls to the ODE function. 
There were 25 calls to the BC function. 
The solution was obtained on a mesh of 100 points.
The maximum residual is  3.046e-06. 
There were 1502 calls to the ODE function. 
There were 26 calls to the BC function. 
The solution was obtained on a mesh of 100 points.
The maximum residual is  5.754e-06. 
There were 1701 calls to the ODE function. 
There were 27 calls to the BC function. 

ans = struct with fields:
    solver: 'bvp4c'
         x: [0 0.1010 0.2020 0.3030 0.4040 0.5051 0.6061 0.7071 0.8081 0.9091 1.0101 1.1111 1.2121 1.3131 1.4141 1.5152 1.6162 1.7172 1.8182 1.9192 2.0202 2.1212 … ] (1×100 double)
         y: [7×100 double]
        yp: [7×100 double]
     stats: [1×1 struct]

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 0.3 0.6 ];
for i =1:numel(rr)
    a= rr(i);
    s=0.5;h=0.1;
    y0 = [1,0,0,0,1,0,0];
    options =bvpset('stats','on','RelTol',1e-4);
    m = linspace(0,10);
    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))];
    figure(2)
    plot(sol.x,(sol.y(6,:)))
    title('Temperature')
    grid on,hold on
    myLegend2{i}=['alfa= ',num2str(rr(i))];
end
figure(1)
legend(myLegend1)
hold on
figure(2)
legend(myLegend2)
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-(2*x+2+1).^2)-(s^2+a1*h^2-(2*x+2+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)))+(((a*2*(x+1)*t)-2*(x+1)-(a*s^2+a*a1*h^2)*t)*du)+(2*a*a3*s*t-a3*s)*dt-(a*a4*s*h+a*a1*s*h)*dt*dv);
        dy(5)=dt;
        dy(6)=ddt;
        dy(7)=(1/(a*(a5*s^2*2*(x+1+1)+a5*h^2*2*(x+1+1))*t-(a5*s^2*2*(x+1+1)+a5*h^2*2*(x+1)-a8*(2*(x+1+1))^3)))*((s^2+h^2+2*a5*h*2-a6*(2*(x+1)).^2-3*a8*2*(x+1+1)*2*(x+1))*ddt);
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]

Thanks alot

[dpb]

@Stephen23 has far more patience (and better eyesight, undoubtedly) than I... <vbg>