Oscillatory solutions problem ODE45

Hi all,

I'm solving this set of ODEs:

function Y=b(zeta,X)

dF1dzeta=X(2);
dF2dzeta=X(1)^2-X(3)^2+X(2)*X(5)+1;
dG1dzeta=X(4);
dG2dzeta=2*X(1)*X(3)+X(4)*X(5);
dH1dzeta=-2*X(1);

Y = [dF1dzeta; dF2dzeta; dG1dzeta; dG2dzeta; dH1dzeta];

Along with the boundary conditions X(1)=X(3)=X(5)=0 at zeta = 0

and X(1)=0, X(3)=1 at zeta = inf

Normally I would use a shooting method such as this

clear all;

dF2=0.0001; 
dG2=0.0001; 

b1=30;
initF2=-0.9420;
initG2=0.7729;

K=zeros(2); 
zetaspan=linspace(0,b1,b1*100+1);
H=[1;1];
options=odeset('AbsTol',1e-8,'RelTol',1e-6);

while max(abs(H))>1e-10,

      [zeta,X]=ode45(@b,zetaspan,[0;initF2+dF2;0;initG2;0],options);
      n=size(zeta,1);
      X1=[X(n,1);X(n,3)-1];
      [zeta,X]=ode45(@b,zetaspan,[0;initF2;0;initG2+dG2;0],options);
      n=size(zeta,1);
      X2=[X(n,1);X(n,3)-1];
      [zeta,X]=ode45(@b,zetaspan,[0;initF2;0;initG2;0],options);
      n=size(zeta,1);
      X3=[X(n,1);X(n,3)-1];

      K(1,1)=(X1(1)-X3(1))/dF2; 
      K(2,1)=(X1(2)-X3(2))/dF2; 
      K(1,2)=(X2(1)-X3(1))/dG2; 
      K(2,2)=(X2(2)-X3(2))/dG2; 

      H=K\-X3;

      initF2=initF2+H(1);
      initG2=initG2+H(2);

end

figure;
hold all;
plot(zeta,X(:,1),'LineWidth',2); 
plot(zeta,X(:,3),'LineWidth',2);  
plot(zeta,X(:,5),'LineWidth',2); 
hold off;
xlabel('$\zeta$','interpreter','latex','FontSize',30);
h=legend('$F$','$G$','$H$','Location','NorthEast');
set(h,'interpreter','latex','FontSize',30)
axis([0 30 -1 2]);
grid on
box on

guessing values for F'(0) and G'(0) to solve the problem. However the solutions are oscillatory as zeta tends to infinity and the above method breaks down for zeta>6 (here zeta=b1). So I reformulated the problem, this time shooting from infinity to zero.

clear all;

dF1=0.00001; 
dG1=0.00001; 
dH1=0.00001; 

b1=30;
initF1=0;
initG1=1;
initH1=1.3494;

K=zeros(3); 
zetaspan=linspace(b1,0,b1*100+1);
H=[1;1;1];
options=odeset('AbsTol',1e-8,'RelTol',1e-6);

while max(abs(H))>1e-10,

      [zeta,X]=ode45(@b,zetaspan,[initF1+dF1;0;initG1;0;initH1],options);
      n=size(zeta,1);
      X1=[X(n,1);X(n,3);X(n,5)];
      [zeta,X]=ode45(@b,zetaspan,[initF1;0;initG1+dG1;0;initH1],options);
      n=size(zeta,1);
      X2=[X(n,1);X(n,3);X(n,5)];
      [zeta,X]=ode45(@b,zetaspan,[initF1;0;initG1;0;initH1+dH1],options);
      n=size(zeta,1);
      X3=[X(n,1);X(n,3);X(n,5)];
      [zeta,X]=ode45(@b,zetaspan,[initF1;0;initG1;0;initH1],options);
      n=size(zeta,1);
      X4=[X(n,1);X(n,3);X(n,5)];

      K(1,1)=(X1(1)-X4(1))/dF1; 
      K(2,1)=(X1(2)-X4(2))/dF1; 
      K(3,1)=(X1(3)-X4(3))/dF1;

      K(1,2)=(X2(1)-X4(1))/dG1; 
      K(2,2)=(X2(2)-X4(2))/dG1;
      K(3,2)=(X2(3)-X4(3))/dG1;

      K(1,3)=(X3(1)-X4(1))/dH1; 
      K(2,3)=(X3(2)-X4(2))/dH1; 
      K(3,3)=(X3(3)-X4(3))/dH1; 

      H=K\-X4;

      initF1=initF1+H(1);
      initG1=initG1+H(2);
      initH1=initH1+H(3);

end

figure;
hold all;
plot(zeta,X(:,1),'LineWidth',2); 
plot(zeta,X(:,3),'LineWidth',2);  
plot(zeta,X(:,5),'LineWidth',2); 
hold off;
xlabel('$\zeta$','interpreter','latex','FontSize',30);
h=legend('$F$','$G$','$H$','Location','NorthEast');
set(h,'interpreter','latex','FontSize',30)
axis([0 30 -1 2]);
grid on
box on

This time using the known values of F(inf)=0 and G(inf)=1 whilst guessing a value for H(inf). But this doesn't seem to get round the problem. What am I best advised to do here? Is there an alternative ODE solver that will combat these oscillatory solutions?