bvp4c failing to produce the correct solution

I'm running the following:

function M = SwirlingFlow(~)  

M = bvpinit(linspace(0,30),@VKinit);

options = bvpset('RelTol',1e-8);

sol = bvp4c(@VK,@VKbc,M,options);

figure;
hold all;
plot(sol.x,sol.y(2,:)); 
plot(sol.x,sol.y(4,:));  
hold off;
xlabel('\zeta')
xlabel('\zeta')
hleg = legend('F\prime','G\prime','Location','NorthEast'); %#ok<NASGU>

figure;
hold all;
plot(sol.x,sol.y(1,:)); 
plot(sol.x,sol.y(3,:)); 
plot(sol.x,(-1)*sol.y(5,:)); 
hold off;
xlabel('\zeta')
hleg = legend('F','G','-H','Location','East'); %#ok<NASGU>

disp('Value of F''(0)')
a = sol.y(2,1);
disp(a)

disp('Value of G''(0)')
b = sol.y(4,1);
disp(b)

disp('Value of H(infinity)')
c = sol.y(5,end);
disp(c)

function yprime = VK(~,y)

yprime = [  y(2) 
            y(1)^2-y(3)^2+y(2)*y(5)
            y(4)
            2*y(1)*y(3)+y(4)*y(5)
            -2*y(1)];

function res = VKbc(ya,yb)

res = [ya(1);ya(3)-1;ya(5);yb(2)-(ya(1)*yb(5));yb(4)-(ya(3)*yb(5))];

function yinit = VKinit(~)

yinit = [0;0.5102;1;-0.6159;0];

I'm sure all the initial values and boundary conditions are correct (I have confirmed the initial conditions using a shooting method). I know exactly how the solution should look and I know that the value of

-H(infty)=0.8845

but using the bvp4c solver it fails to produce the correct solution. Can anybody help debug this?

Thanks.