why my theta is equal to zero?

tryanderror123()

ans = 1×2

8 5

First solution:
-f"(0) = 1.048998010
f`(0) = 0.580400796
F`(0) = 0.268610539
-theta`(0) = -0.000000000
-theta_p(0) = -0.000000000

function tryanderror123
format long g 
%Define all parameters 
global Pr beta rho_r gamma A C SS M k 
Pr = 7; %Prandtl number
beta = 0.5; %fluid particle interaction parameter
rho_r = 10; %relative density particle and fluid
gamma = 0.8; %ratio of specific heat of the fluid and dust phase
A = 0.4; %slip parameter
C = 1; % stretching/shrinking
SS = 0; %suction parameter
M = 0; %magnetic parameter
k = 0; %porosity parameter
%Boundary layer thickness & stepsize
etaMin = 0;
etaMax1 = 5;
etaMax2 = 5; %15, 10
stepsize1 = etaMax1;
stepsize2= etaMax2;
%%%%%%%%%%%%%%%%%%%%%%   first solution   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
options = bvpset('stats','off','RelTol',1e-10); 
solinit = bvpinit (linspace (etaMin, etaMax1, stepsize1), @OdeInit1);
sol = bvp4c(@OdeBVP, @OdeBC, solinit, options);
eta = linspace (etaMin, etaMax1, stepsize1);
y = deval (sol, eta);
%Displaying the output for first solution 
fprintf('\nFirst solution:\n');
fprintf('-f"(0) = %7.9f\n',-y(3));
fprintf('f`(0) = %7.9f\n', y(2))
fprintf('F`(0) = %7.9f\n', y(7))
fprintf('-theta`(0) = %7.9f\n',-y(5));
fprintf('-theta_p(0) = %7.9f\n',-y(8));
fprintf('\n');
end
%Define the ODE function
function ff = OdeBVP(x, y, Pr, beta, rho_r, gamma, M, k)
global Pr beta rho_r gamma M k
ff = [
    y(2) %f'
    y(3) %f''
    y(2)*y(2)-y(1)*y(3)+M*y(2)+k*y(2)-beta*rho_r*(y(7)-y(2)) %f'''
    y(5) %theta'
    Pr*( 2*y(2)*y(4) - y(1)*y(5) - ((2*beta*rho_r)/(3*Pr))*(y(8)-y(4)) ) %theta''
    y(7) %F'
    ( y(7)*y(7) - beta*(y(2)-y(7)) ) / y(6) %F''
    ( 2*y(7)*y(8) + ((2*gamma*beta)/(3*Pr))*(y(8)-y(4)) ) / y(6) %theta_p'
    ]; 
end
%Define the boundary condition
function res = OdeBC(ya, yb, A, C, SS) 
global A C SS
res = [
    ya(2)-C-A*ya(3)
    ya(1)-SS
    ya(4)
    yb(2)
    yb(7)
    yb(6)-yb(1)
    yb(4)
    yb(8) ];
end
%Setting the initial guess for first solution  
function v = OdeInit1 (x, A, C, SS)
global A C SS 
v = [0
    0
    0
    0
    0
    0.3
    0
    0];  
end

Answer 1

[Walter Roberson]

deval() returns a matrix with one row for each state variable, and one column for each time evaluated at.

You are using linear indexing with small indices, so you are extracting data from the first column -- which is the column associated with the initial time. So your first column should correspond to the initial conditions that were used. It should not be surprsing if those come out as 0.