it means the size of y variable is 8 abd you are trying to get 9th 10th 11th index of the y vector. This is because your initial solution, defined in sol_init variable consist of only 8 elements, where it should be 11. So a simple fix is you fix this line:
%solinit = bvpinit(linspace(0, 2, 50), [1 0 0 0 1 1 1 1]);
Then you had some syntax error in bcexam2() function. When you put space between x +T in an array element, matlab can think it is a 2 element vector. You should use paranthesis or no space when you define that way.
And lastly, since your solution is 11x1 vector, you should modify your bcexam2 function to provide 11x1 vector instead of 9x1. I cannot solve this one as I have no idea what equation you are trying to solve.
Latest_code
function Latest_code
%%
for Pr = [1 2 3 4]
% Define constants
M = 0.5;
B = 0.5;
C = 0.2;
A = 0.1;
R = 0.1;
E = 0.1;
K = 0.5;
W = 0.1;
N = 0.1;
T = 0.3;
Q = 0.1;
N1 = 0.1;
N2 = 0.1;
S = 0.1;
d = 1;
m = 1;
E1 = 0.1;
Pe = 0.2;
S1 = 0.1;
M1 = 0.1;
O = 0.1;
% Defining initial guess
solinit = bvpinit(linspace(0, 2, 50), [1 0 0 0 1 1 1 1 1 1 1]);
% Solving the boundary value problem
sol = bvp4c(@bvpexam2, @bcexam2, solinit);
x = sol.x;
y = sol.y;
% Plotting of velocity
figure(2)
plot(x, y(4, :), 'linewidth', 1)
ylim([0 1])
xlim([0 9])
hold on
xlabel('\eta', 'FontWeight', 'bold', 'FontSize', 16)
ylabel('\theta(\eta)', 'FontWeight', 'bold', 'FontSize', 16)
end
% System of first order ODEs
function yvector = bvpexam2(~, y)
yy1 = (1 / (1 + (1 / B))) * (2 * y(2) * y(2) - y(1) * y(3) + M * y(2) + K * y(2) - W * (y(4) - N * y(7) - R));
yy2 = -(1 / (((4 / 3) * R) + 1)) * (Pr * (y(1) * y(5) - y(2) * y(4) - T * y(2) + Q * y(4) + N1 * y(8) * y(5) + N2 * y(5) * y(5) + E * (1 + 1 / B) * y(3) * y(3) + (E * y(2) * y(2)) * (M + K)));
yy3 = -S * (y(1) * y(8) - y(2) * y(7) - C * y(2) - A * y(7) * ((1 + d * y(4)) ^ m) * exp(-E1 / (1 + d * y(4)))) - (N2 / N1) * y(6);
yy4 = (Pe * (y(9) * (y(10) + O) + y(11) * y(8))) - (S1 * ((y(11) - y(10)) * y(1))) + (S1 * M1 * y(2));
yvector = [ y(2); y(3); yy1; y(5); yy2; y(8); yy3 ; y(11); yy4 ];
end
% Residual of the boundary conditions
function residual = bcexam2(y0, yinf)
residual = [(y0(1)- 1); (y0(2) - 1); (y0(4) - 1 +T) ; (y0(7)-1+C) ; (y0(10)-1+M1) ; yinf(2); yinf(4); yinf(7); yinf(10)];
end
end
