Please help me to run this simple code

% Error

Array indices must be positive integers or logical values.

Error in Untitled2 (line 45)

kh=kf*((k1+2*kf-2*ph1(kf-k1))/(k1+2*kf+ph1*(kf-k1)));

% code
proj()

Array indices must be positive integers or logical values.

Error in solution>proj (line 47)
kh=kf*((k1+2*kf-2*ph1(kf-k1))/(k1+2*kf+ph1*(kf-k1)));
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

function sol = proj
    clc; clf; clear;
    rhof=997.1*10^-3;kf=0.613*10^5;cpf=4179*10^4;muf=10^-3*10;
alfaf=kf/(rhof*cpf);

bef=21*10^-5;
sigf=0.05*10^-8;
ky=muf/rhof;
%Titanium oxide
ph4=0.01;
rho4=4250*10^-3;
cp4=711*10^4;
k4=8.953*10^5;

sig4=2.6*10*10^-1;

%Ag
ph1=0.01;
rho1=10500*10^-3; 
cp1=235*10^4;
k1=429*10^5;be1=21*10^-5;
sig1=0.74*10^-2;

%copper
ph2=0.01;
rho2=8933*10^-3;
cp2=385*10^4;
k2=400*10^5;
sig2=5.96*10^-1;
be2=1.67*10^-5;

%Alumina
ph3=0.01;
rho3=3970*10^-3;
cp3=765*10^4;
k3=40*10^5;
be3=0.85*10^-5;
sig3=3.5*10^-1;


%Relation of ternary hyprid


kh=kf*((k1+2*kf-2*ph1(kf-k1))/(k1+2*kf+ph1*(kf-k1)));


kh=kn*((k2+2*kn-2*ph2(kn-k2))/(k2+2*kn+ph2*(kn-k2)));

kt=kh*((k3+2*kh-2*ph3(kh-k3))/(k3+2*kh+ph3*(kh-k3)));

kT=kt*((k4+2*kt-2*ph4*(kt-k4))/(k4+2*kt+ph4*(kt-k4)));

muT= muf/((1-ph1)^2.5*(1-ph2)^2.5*(1-ph3)^2.5*(1-ph4)^2.5);

rhoT=rhof*((1-ph4)*((1-ph3)*((1-ph2)*((1-ph1)+ph1*((ph1*rho1*(1/rhof))))+(ph2*rho2*(1/rhof))))+(ph3*rho3*(1/rhof))+ph4*rho4*(1/rhof));

sign = sigf*(1+(3*(sigf-1)*ph1)/((sigf+2)-(sigf-1)*ph1));
sigh = sign*(1+(3*(sign-1)*ph2)/((sign+2)-(sign-1)*ph2));
sigt = sigh*(1+(3*(sigh-1)*ph3)/((sigh+2)-(sigh-1)*ph3));
sigT = sigt*(1+(3*(sigt-1)*ph4)/((sigt+2)-(sigt-1)*ph4));


%vt=rhot*cpt
VT=(rhof*cpf)*((1-ph4)*((1-ph3)*((1-ph2)*((1-ph1)+ph1*((ph1*cp1*(1/rhof*cpf))))+(ph2*cp2*(1/rhof*cpf))))+(ph3*cp3*(1/rhof*cpf))+ph4*cp4*(1/rhof*cpf));

%disp('vt');disp(vt);
%vb=rho*betb
vT =(1-ph1)*((1-ph2)*((1-ph3)+ph3*((rho3*be3)/(rhof*bef)))+ph2*((rho2*be2)/(rhof*bef)))+(1-ph1)*ph1*((rho1*be1)/(rhof*bef));  
 myLegend1 = {};myLegend2 = {};
rr = [0.1 0.3 0.5 0.7]
for i =1:numel(rr)
    Rd= rr(i)
M=0.5;
    R=1;Pr=6.9;
    m = linspace(0,1);
y0 = [0,1,0,1,0,0,0,0];options =bvpset('stats','on','RelTol',1e-5);
   solinit = bvpinit(m,y0);
      sol= bvp4c(@projfun,@projbc,solinit,options);
          figure(1)
    plot(sol.x,(sol.y(1,:)))
    grid on,hold on
    myLegend1{i}=['n= ',num2str(rr(i))];
    
    figure(2)
 plot(sol.x,(sol.y(7,:)))
    grid on,hold on
    myLegend2{i}=['n= ',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(8,1);
        % alignComments
        f = y(1);
        df = y(2);
         g = y(3);
        dg= y(4);
        h= y(5);
        dh = y(6);
         t = y(7);
       dt=y(8);
               dy(1) = df; 
        dy(2) =(1/(1+x^2))*(3*x*df+R*((muT/muf))*(f^2-f*df+h*df-g^2+(sigT/sigf)*M*f));
         dy(3) = dg;
     dy(4) = (1/(1+x^2))*(3*x*dg+R*((muT/muf))*(2*f*g-f*dg+h*dg+(sigT/sigf)*M*g));
         dy(5) =dh ;
      dy(6) =(1/(1+x^2))*(2*x*dh-h+R*((muT/muf))*(f*dh-h*dh-f*h));
         dy(7) =dt;
       dy(8)=(1/(1+x^2+(4/3)*Rd*(1/(kT/kf))))*(((x-4)*dt)-4*t-Pr*R*((vT/(rhof*cpf))*(kf/kT))*(x*f*dt-2*f*t-h*dt));
                     

    end
end

function res= projbc(ya,yb)
res= [ya(1)-0.1;
  ya(3)-1;
  ya(5)-0.1;
    ya(7)-1;
    
    yb(1)-0.1;
        yb(3);
        yb(5);
     yb(7);
       ];
       end

[Torsten]

Use

kh=kf*((k1+2*kf-2*ph1*(kf-k1))/(k1+2*kf+ph1*(kf-k1)));

instead of

kh=kf*((k1+2*kf-2*ph1(kf-k1))/(k1+2*kf+ph1*(kf-k1)));

[T K]

I am uses this above relation

kh=kf*((k1+2*kf-2*ph1*(kf-k1))/(k1+2*kf+ph1*(kf-k1)));

but the code not Run??

[T K]

Please help me, Prof Torsten

[Torsten]

Did you do what I suggested ?

Your line misses a multiplication sign - I added it.

[T K]

Iam added multiplication sign but the code not run.Pleae help me

% Error using bvp4c (line 251)

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

Error in proj (line 73)

sol= bvp4c(@projfun,@projbc,solinit,options);

% code

proj()

rr = 1×4

0.1000 0.3000 0.5000 0.7000

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

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

Rd = 0.1000

Warning: Unable to meet the tolerance without using more than 1250 mesh points.
The last mesh of 892 points and the solution are available in the output argument.
The maximum residual is 0.240883, while requested accuracy is 1e-05.

Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 9.13308e-94.

The solution was obtained on a mesh of 2674 points.
The maximum residual is  2.409e-01. 
There were 107407 calls to the ODE function. 
There were 354 calls to the BC function. 

Rd = 0.3000

Warning: Unable to meet the tolerance without using more than 1250 mesh points.
The last mesh of 892 points and the solution are available in the output argument.
The maximum residual is 0.240883, while requested accuracy is 1e-05.

Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 1.43969e-92.

The solution was obtained on a mesh of 2674 points.
The maximum residual is  2.409e-01. 
There were 107374 calls to the ODE function. 
There were 354 calls to the BC function. 

Rd = 0.5000

Warning: Unable to meet the tolerance without using more than 1250 mesh points.
The last mesh of 892 points and the solution are available in the output argument.
The maximum residual is 0.240883, while requested accuracy is 1e-05.

Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 1.26505e-91.

The solution was obtained on a mesh of 2674 points.
The maximum residual is  2.409e-01. 
There were 107369 calls to the ODE function. 
There were 354 calls to the BC function. 

Rd = 0.7000

Warning: Unable to meet the tolerance without using more than 1250 mesh points.
The last mesh of 892 points and the solution are available in the output argument.
The maximum residual is 0.240883, while requested accuracy is 1e-05.

Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 9.33478e-91.

The solution was obtained on a mesh of 2674 points.
The maximum residual is  2.409e-01. 
There were 107331 calls to the ODE function. 
There were 354 calls to the BC function. 

ans = struct with fields:
    solver: 'bvp4c'
         x: [0 0.0011 0.0022 0.0034 0.0045 0.0056 0.0067 0.0079 0.0090 0.0101 0.0112 0.0123 0.0135 0.0146 0.0157 0.0168 0.0180 0.0191 0.0202 0.0213 0.0224 0.0236 … ] (1×892 double)
         y: [8×892 double]
        yp: [8×892 double]
     stats: [1×1 struct]

function sol = proj
    clc; clf; clear;
    rhof=997.1*10^-3;kf=0.613*10^5;cpf=4179*10^4;muf=10^-3*10;
alfaf=kf/(rhof*cpf);

bef=21*10^-5;
sigf=0.05*10^-8;
ky=muf/rhof;
%Titanium oxide
ph4=0.01;
rho4=4250*10^-3;
cp4=711*10^4;
k4=8.953*10^5;

sig4=2.6*10*10^-1;

%Ag
ph1=0.01;
rho1=10500*10^-3; 
cp1=235*10^4;
k1=429*10^5;
be1=21*10^-5;
sig1=0.74*10^-2;

%copper
ph2=0.01;
rho2=8933*10^-3;
cp2=385*10^4;
k2=400*10^5;
sig2=5.96*10^-1;
be2=1.67*10^-5;

%Alumina
ph3=0.01;
rho3=3970*10^-3;
cp3=765*10^4;
k3=40*10^5;
be3=0.85*10^-5;
sig3=3.5*10^-1;


%Relation of ternary hyprid
kn=kf*((k1+2*kf-2*ph1*(kf-k1))/(k1+2*kf+ph1*(kf-k1)));
kh=kn*((k2+2*kn-2*ph2*(kn-k2))/(k2+2*kn+ph2*(kn-k2)));

kt=kh*((k3+2*kh-2*ph3*(kh-k3))/(k3+2*kh+ph3*(kh-k3)));

Km=kt*((k4+2*kt-2*ph4*(kt-k4))/(k4+2*kt+ph4*(kt-k4)));

muT= muf/((1-ph1)^2.5*(1-ph2)^2.5*(1-ph3)^2.5*(1-ph4)^2.5);

rhoT=rhof*((1-ph4)*((1-ph3)*((1-ph2)*((1-ph1)+ph1*((ph1*rho1*(1/rhof))))+(ph2*rho2*(1/rhof))))+(ph3*rho3*(1/rhof))+ph4*rho4*(1/rhof));

sign = sigf*(1+(3*(sigf-1)*ph1)/((sigf+2)-(sigf-1)*ph1));
sigh = sign*(1+(3*(sign-1)*ph2)/((sign+2)-(sign-1)*ph2));
sigt = sigh*(1+(3*(sigh-1)*ph3)/((sigh+2)-(sigh-1)*ph3));
sigT = sigt*(1+(3*(sigt-1)*ph4)/((sigt+2)-(sigt-1)*ph4));


%vt=rhot*cpt
vt=(rhof*cpf)*((1-ph4)*((1-ph3)*((1-ph2)*((1-ph1)+ph1*((ph1*cp1*(1/rhof*cpf))))+(ph2*cp2*(1/rhof*cpf))))+(ph3*cp3*(1/rhof*cpf))+ph4*cp4*(1/rhof*cpf));

 myLegend1 = {};myLegend2 = {};
rr = [0.1 0.3 0.5 0.7]
for i =1:numel(rr)
    Rd= rr(i)
M=0.5;
    R=1;Pr=6.9;
    m = linspace(0,1);
y0 = [0,1,0,1,0,0,0,0];options =bvpset('stats','on','RelTol',1e-5);
   solinit = bvpinit(m,y0);
      sol= bvp4c(@projfun,@projbc,solinit,options);
          figure(1)
    plot(sol.x,(sol.y(1,:)))
    grid on,hold on
    myLegend1{i}=['n= ',num2str(rr(i))];
    
    figure(2)
 plot(sol.x,(sol.y(7,:)))
    grid on,hold on
    myLegend2{i}=['n= ',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(8,1);
        % alignComments
        f = y(1);
        df = y(2);
         g = y(3);
        dg= y(4);
        h= y(5);
        dh = y(6);
         t = y(7);
       dt=y(8);
               dy(1) = df; 
        dy(2) =(1/(1+x^2))*(3*x*df+R*((muT/muf))*(f^2-f*df+h*df-g^2+(sigT/sigf)*M*f));
         dy(3) = dg;
     dy(4) = (1/(1+x^2))*(3*x*dg+R*((muT/muf))*(2*f*g-f*dg+h*dg+(sigT/sigf)*M*g));
         dy(5) =dh ;
      dy(6) =(1/(1+x^2))*(2*x*dh-h+R*((muT/muf))*(f*dh-h*dh-f*h));
         dy(7) =dt;
       dy(8)=(1/(1+x^2+(4/3)*Rd*(1/(Km/kf))))*(((x-4)*dt)-4*t-Pr*R*((vt/(rhof*cpf))*(kf/Km))*(x*f*dt-2*f*t-h*dt));         

    end
end

function res= projbc(ya,yb)
res= [ya(1)-0.1;
  ya(3)-1;
  ya(5)-0.1;
    ya(7)-1;
    
    yb(1)-0.1;
        yb(3);
        yb(5);
     yb(7);
       ];
       end

[Torsten]

The warning message of the recent MATLAB release is different (see above).

MATLAB cannot solve your equations. Maybe there are mistakes in the expressions or in the parameters, maybe the equations are just very difficult to solve - we cannot tell.