Numerical Solution for Nonlinear Shooting Method

I am having issues with the code for the numerical solution of the boundary problem. I have an error somewhere in my code that is shifting the graph of the numerical solution up an away from the exact solution. It is something I am overlooking, can somebody help locating the error?

% Nonlinear Shooting Method Example using Euler method
  % Inputs: interval inter, initial vector y0, number of steps n
  % Output: time steps t, solution y
  % Example usage: NLShooting([1 3],[17 43/3],20);
  function NLShooting(inter,bc,n)
  alpha=bc(1);
  beta=bc(2);
  MaxSecantIterations=60;
  F=zeros(1,MaxSecantIterations);
  t=zeros(1,n);
  y=zeros(n,2);
  t(1)=inter(1);
  h=(inter(2)-inter(1))/n;
  y(1,1)=alpha;
  s(1)=(beta-alpha)/(inter(2)-inter(1));
  y(1,2)=s(1);
  for i=1:n
    t(i+1)=t(i)+h;
    y(i+1,:)=eulerstep(t(i),y(i,:),h);
  end
  F(1)=y(n,1)-beta;
  s(2)=s(1)+(F(1))/(inter(2)-inter(1));
  y(1,2)=s(2);
  j=2;
  while (abs(F(j-1))>1/1000) && j<MaxSecantIterations
  for i=1:n
    t(i+1)=t(i)+h;
    y(i+1,:)=eulerstep(t(i),y(i,:),h);
  end
  F(j)=y(n,1)-beta;
  s(j+1)=s(j)-F(j)*(s(j)-s(j-1))/(F(j)-F(j-1));
  y(1,2)=s(j+1);
  y(1,1)=alpha;
  j=j+1;
  y(:,1)
  end
  plot(t,y(:,1),'go');
  hold on
  plot(t,t.^2+16./t,'blue');
  grid
  title('Nonlinear Shooting Method');
  function y=eulerstep(t,y,h)
  %one step of the Euler method
  %Input: current time t, current vector y, step size h
  %Output: the approximate solution vector at time t+h
  y=y+h*ydot(t,y);

    function z=ydot(t,y)
    z(1) = y(2);
    z(2) =  1/8*(32+2*t^3-y(1)*y(2));