How to implement the shooting method

Hi. I'm attempting to solve the the ODEs one arrives at when analysing the problem of the flow due to a rotating disk, that is:

F''=F^2+G^2+F'H
G''=2FG+G'H
H''=-2F'

with BCs

F(0)=0, 
G(0)=1, 
H(0)=0 
F(infinity)=G(infinity)=0, 
H(infinity)=-A (constant)

I have implemented the BVP solver as such:

function M = Script(~) 

M = bvpinit(linspace(0,4.5,801),@VKinit);
sol = bvp4c(@VK,@VKbc,M);

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)
            y(6)
            -2*y(2)];

% y(1) = F, y(2) = F', y(3)= G, y(4) = G', y(5) = H and y(6)= H'.

function res = VKbc(ya,yb)

res = [ya(1);ya(3)-1;ya(5);yb(1);yb(3);yb(5)+0.88447];

function yinit = VKinit(~)

yinit = [0;0;0;0;0;0];

Plotting the solutions to this gives me exactly what I want, but that's only because I have taken the correct value of H(infinity)=-0.88447 from previous texts. (In this case I have taken infinity=4.5, the solution converges by this point)

I need to be able to find the value of this unknown boundary condition myself.

How do I go about modifying my code to obtain this value of H(infinty) for myself using a shooting method? Is the shooting method in fact the best way of evaluating H(infinity)?

Any help and advice anyone could give would be greatly appreciated!