I need curve for f "(0) vs phi with G= 0; 0.1; 0.2; 0.3

function main

p1 = 0:0.02:2;

for k = 1:length(p1)

phi = p1(k);

end

Pr=6.2;

G=input('G='); % G=0; 0.1; 0.2; 0.3

rhos=997.1;Cps=4179;ks=0.613; rhof=8933;Cpf=385;kf=401;

a1=(1-phi)^2.5*(1-phi+phi*(rhos/rhof));a2=(1-phi+phi*((rhos*Cps)/(rhof*Cpf)));A=(ks+2*kf+phi*(kf-ks))/(ks+2*kf-2*phi*(kf-ks));

xa=0;xb=5;

solinit=bvpinit(linspace(xa,xb,100),[1 1 0 1 1]);

sol=bvp4c(@ode,@bc,solinit);

xint=linspace(xa,xb,101);

sxint=deval(sol,xint);

function res=bc(ya,yb)

res=[ya(1); ya(2)-1-G*ya(3); ya(4)-1; yb(2); yb(4)];

end

function dydx=ode(x,y)

dydx=[y(2); y(3); a1*(y(2)^2-y(3)*y(1)); y(5); -A*Pr*a2*y(1)*y(5)];

end

f0 = deval(sol,0);

f0(3);

figure(1)

plot(p1,sxint([3],:));

hold on

end