With the given program, I need to draw B & C (defined in the program) with PHI variations (0, 0.1, 0.2) ANYBODY HELP and write the code to draw the required fig thanks in Advance

MATLAB PROGRAM

%I need to draw B & C (defined in the program) with PHI variations (0, 0.1, 0.2)
function main
format('long');
gg=['r','k','b','g','m','c','y','r.','m.','k.'];
A=0.5;  Pr=7;   p=0.5;    
% phi=0.1;
phi=input('phi=');
rhos=997;Cps=4179;ks=0.613;             %for    WATER
rhof=385;Cpf=8933;kf=400;               %for    Cu
a1=(1-phi)^2.5*(1-phi+phi*(rhos/rhof));
a2=(1-phi)^2.5*(1-phi+phi*((rhos*Cps)/(rhof*Cpf)));
Knf=(kf)*(ks+2*kf-2*phi*(kf-ks))/(ks+2*kf+phi*(kf-ks));

% B=-(Knf/kf)*T'(0);            f=y(1),f'=y(2),f"=y(3),T=y(4),T'=y(5)
% C=f"(0)/(1-phi)^2.5;
%%%%%%%%%%%%%%%%%%%

xa=0;xb=6;
solinit=bvpinit(linspace(xa,xb,1000),[0 0 0 1 1 0 1 0]);
sol=bvp4c(@ode,@bc,solinit);
xint=linspace(xa,xb,100);
sxint=deval(sol,xint);

%Boundary Condition
function res=bc(ya,yb)
res=[ya(1); ya(2)-1; ya(4); ya(5)-p; ya(7)-1; yb(2); yb(5); yb(7)];        
end 

function dydx=ode(x,y)    
dydx=[y(2);  y(3);  2*a1*y(2)*(y(2)+y(5))-a1*y(3)*(y(1)+y(4));  
      y(5);  y(6);  2*a1*y(5)*(y(2)+y(5))-a1*y(6)*(y(1)+y(4));
      y(8);  A*Pr*a2*y(7)*(y(2)+y(5))-Pr*a2*y(8)*(y(1)+y(4))];              
end

plot(xint,sxint([2],:));    %for f'
xlabel('\eta'); 
ylabel('f`');
hold on

[Walter Roberson]

Note: your gg should be

gg = {'r','k','b','g','m','c','y','r.','m.','k.'};

However you never use gg in the code you posted, so you would not have noticed.

[MINATI PATRA]

Hi Walter I saw late Thanks for your response. Actually I need the plot for 'B' & 'C'(defined in the program) vs phi(ranges from 0 to 0.2) WHERE already 'x' Varies previously. I have made Other figs

[MOSLI KARIM]

Here is the correction to your code

function ANSWER5
format long
global A Pr p
col={'r','k','b'};
A=0.5;  Pr=7;   p=0.5;    
phi=[0,0.1,0.2 ];          

for ik=1:length(phi)
    rhos=997;Cps=4179;ks=0.613;             %for    WATER
    rhof=385;Cpf=8933;kf=400;               %for    Cu
    a1=(1-phi(ik))^2.5*(1-phi(ik)+phi(ik)*(rhos/rhof));
    a2=(1-phi(ik))^2.5*(1-phi(ik)+phi(ik)*((rhos*Cps)/(rhof*Cpf)));
    Knf=(kf)*(ks+2*kf-2*phi(ik)*(kf-ks))/(ks+2*kf+phi(ik)*(kf-ks));
    
    % B=-(Knf/kf)*T'(0);            f=y(1),f'=y(2),f"=y(3),T=y(4),T'=y(5)
    % C=f"(0)/(1-phi)^2.5;
    %%%%%%%%%%%%%%%%%%%
    
    xa=0;xb=6;
    solinit=bvpinit(linspace(xa,xb,1000),[0 0 0 1 1 0 1 0]);
    sol=bvp4c(@ode,@bc,solinit);
    xint=linspace(xa,xb,100);
    sxint=deval(sol,xint);
    
    figure(1)
    hold on 
    plot(xint,sxint(2,:),col{ik} );    %for f'
    xlabel('\eta'); 
    legend('\phi=0','\phi=0.1','\phi=0.2')
end 
%Boundary Condition
    function res=bc(ya,yb)
        res=[ya(1); 
            ya(2)-1; 
            ya(4); 
            ya(5)-p; 
            ya(7)-1; 
            yb(2);
            yb(5); 
            yb(7)];        
    end 

    function dydx=ode(x,y)    
        dydx=[y(2);  
            y(3); 
            2*a1*y(2)*(y(2)+y(5))-a1*y(3)*(y(1)+y(4));  
            y(5); 
            y(6);  
            2*a1*y(5)*(y(2)+y(5))-a1*y(6)*(y(1)+y(4));
            y(8);  
            A*Pr*a2*y(7)*(y(2)+y(5))-Pr*a2*y(8)*(y(1)+y(4))];              
    end


end