clear all
close all

% where we solve
a = 
b = 

% parameters for the Newton's method
eps = 1e-6;
maxiter = 100;
err = 1;
k = 1;

% initial guess
eta = 

% repeat until the second BC is not satisfied
while (err >= eps && k < maxiter)

    % with the given IC solve the 4 ODE (using R-K methods)
    ic = [y10, y20, p10, p20];
    [x, y] = ode45(@varode, [a, b], ic);    

    % compute the new aproximation of eta by the Newton's method
    eta_new = eta - 
    
    err = abs(eta_new - eta); % error
    eta = eta_new; % use the new eta in the next step
    k = k+1; 

    fprintf('k = %d, eta = %10.8f, error = %10.8e\n', k-1, eta, err);
end

if (err > eps)
    fprintf("Shooting method does not converge.\n");
    return;
end

figure
plot (x, y(:,1),'Linewidth',2);
legend('y1');
xlabel('x');
ylabel('y');
figure
plot (x, y(:,2),'Linewidth',2);
legend('y2');
xlabel('x');
ylabel('y');

% -----------------------------------------------------------------------------------
    
function dy = varode(x,y)

    dy = zeros(size(y));

    % 4 ODE - original ODE and variational ODE
    % the solution has 4 components:
    % y(1) - y1, y(2) - y2, y(3) - p1, y(4) - p2

    dy(1) = 
    dy(2) = 
    dy(3) = 
    dy(4) = 
end