format compact
format long

tol = 0.000000000001;
a = 1;
b = 3;
xspan = [a b];
alpha = 17;
beta = 43/3;
t = (beta - alpha)/(b - a);
plot(b, beta, 'r*');
hold on
ylim([10, 21]);
title('Nonlinear shooting method');
xlabel('x');
ylabel('y');
grid on

k = 0;
while true

    y0 = [alpha, t, 0, 1];
    opts = odeset('RelTol',1e-10,'AbsTol',1e-10);
    [x, y] = ode45(@nleq, xspan, y0);

    k = k + 1
    error = abs(y(end, 1) - beta)
    t = t - (y(end, 1) - beta)/y(end, 3);

    if error < tol
        break
    end

    plot(x, y(:, 1))

    pause(3);
end

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

