Return them as function outputs:
[x,y] = Shrinking_bvp4c % here are the values
plot(x, y(2, :) ,'linewidth', 1)
hold on
xlabel('\eta', 'fontweight', 'bold', 'fontsize', 16)
ylabel('f^/(\eta)', 'fontweight', 'bold', 'fontsize', 16)
function [x,y] = Shrinking_bvp4c % <- define the output arguments here.
%!!!!!!!!!!! Do NOT put anti-pattern CLEAR at the top of a function !!!!!!!!!!
% defining parameters
k=0.2;
K=0.2;
M=0.2;
S=2;
sol1 = bvpinit(linspace(0,3,30),[1 0 1 0]);
sol = bvp4c(@bvp2D,@bc2D,sol1);
x = sol.x;
y = sol.y;
%% Residual of the boundary conditions
function residual = bc2D(y0, yinf)
residual=[y0(1)-S; y0(2)+1; yinf(2); yinf(3)];
end
%% System of First Order ODEs
function yvector = bvp2D(t,y)
yy1 = 1/y(1)*(y(2)*y(2)-y(1)*y(3)-y(4)+k*(2*y(2)*y(4)-y(3)*y(3))+M*y(2)+K*y(2));
yvector = [y(2);y(3);y(4);yy1];
end
end
