Ths usual solution for that is to simply increase the initial mesh size, and here that is:
z = linspace(0,1,7500);
It then runs without error. (I used yyaxis because 
is more than an order of magnitude greater than the others, and those details get lost plotting them on the same axes.)
is more than an order of magnitude greater than the others, and those details get lost plotting them on the same axes.) Also, please never use global variables! Pass them as extra parameters, as I did here.
Increasing the size of ‘z’, correcting the global variable problem, and plotting with yyaxis are the only changes I made to the posted code.
% global Pem Peh Da b U yw yf
Pem = 2000; Peh = 500;
Da = 1e8; b = 0.015;
U = 1; yw = 0.9; yf = 0.0503;
z = linspace(0,1,7500);
%guess = 4*ones(1,2);
solinit = bvpinit(z,[1,0,1,0]);
sol = bvp4c(@(t,y)myfun(t, y, Pem,Peh, Da, b, U, yw, yf),@(ya,yb)mybc(ya,yb, Pem, Peh, yf),solinit)
figure
yyaxis left
plot(sol.x, sol.y(1:3,:))
yyaxis right
plot(sol.x, sol.y(4,:))
grid
legend(compose('y_%d',1:size(sol.y,1)), 'Location','S')
function F = myfun(t, y, Pem,Peh, Da, b, U, yw, yf)
% global Pem Peh Da b U yw
y1 = y(1);
y2 = y(2);
y3 = y(3);
y4 = y(4);
F(1) = y2;
F(2) = Peh*(y2 + b*Da*y3*exp(-1/y1) - U*(yw-y1));
F(3) = y4;
F(4) = Pem*(y4 + Da*y3*exp(-1/y1));
F = F(:);
end
function f = mybc(ya,yb, Pem, Peh, yf)
% global Pem Peh yf
f(1) = ya(2)+Peh*(yf-ya(1));
f(2) = yb(2);
f(3) = ya(4) + Pem*(1-ya(3));
f(4) = yb(4);
f=f(:);
end
function g = guess(z) % initial guess for y and y'
g = [0.5
0.5
0.5
0.5];
end
Experiment to get different results.
.
