Solving ODEs using bvp4c

I am looking to solve the following ODEs with the given initial and boundary conditions.(ref: Journal of Engineering Science and Technology Vol. 5, No. 2 (2010) 190 - 212; NUMERICAL SOLUTION OF STEADY STATE DISPERSION FLOW MODEL FOR LACTOSE-LACTASE HYDROLYSIS WITH DIFFERENT KINETICS IN FIXED BED). The method used in the mentioned paper is 'Orthogonal Collocation'.

I do not know how to input the sqrt(u) term in the boundary condition.

I have attached the piece of code and would appreciate your assistance. Thank you.

%% Main program
clc,clear;close all;
global Pe 
Pe = 5; % Peclet number
xmesh = linspace(0,1); % discretisation of the 1D domain
solinit = bvpinit(xmesh,[0;0]); % initial guess

sol = bvp4c(@biofcn,@biobc,solinit); % calling the inbuilt solver

% Plotting the results
plot(sol.x,sol.y(1,:))
title('substrate concentration in a tubular bioreactor')
xlabel('dimensionless axial coordinate')
ylabel('dimensionless concentration')

%% functions
function dydx = biofcn(x,y)
global Pe
phi = 0.74; phithree = 3.42;Cs0 = 0.292;Km = 0.00156;
dydx = [y(2)
        (phithree*(Cs0*y(1)+phi+Km-sqrt((phi+Km-Cs0*y(1)).^2+4*Km*Cs0*y(1)))+2*sqrt(x)*y(2))*Pe/4*x];
end

function res = biobc(ya,yb)
global Pe
res = [(2*ya(2)/Pe)+1-ya(1);yb(2)];
end