tol = 0.000001;
xspan = [1 3];
alpha = 17;
beta = 43/3;
T = (beta - alpha)/(xspan(2) - xspan(1));
hold on
ylim([10 22]);
title('Non-linear shooting method');
xlabel('x');
ylabel('y');

while true
    y0 = [alpha T 0 1];

    [x, y] = ode45(@nleq, xspan, y0);
    T = T - (y(end, 1) - beta)/y(end, 3);
    if (abs(y(end, 1) - beta) < tol)
        break;
    end
    
    T
    abs(y(end, 1) - beta)

    plot(x, y(:, 1));
    pause(3);

end

function dydx = nleq(x, y)
%NLEQ Evaluate the rhs of the ODE
    dydx(1, 1) = y(2);
    dydx(2, 1) = -y(1)*y(2)/8 + x^3/4 + 4;
    dydx(3, 1) = y(4); 
    dydx(4, 1) = -y(3)*y(2)/8 - y(4)*y(1)/8;
end

