function [x, theta2] = solutions_Final_p2b(theta1, L1, L2, h)

    x = zeros(size(theta1));
    theta2 = x;
    
    for i=1:length(theta1)
        [x(i), theta2(i)] = solutions_Final_p2b_Newton(theta1(i), L1, L2, h);
    end
    
end

function [x, theta2] = solutions_Final_p2b_Newton(theta1, L1, L2, h)
    
    % Absolute tolerance
    tol = 1E-8;

    % Initial guess    
    x = 1;
    theta2 = 0.0;
    sol = [x;theta2];
    
    % Initial residual
    r = solutions_Final_p2b_Residual(x, theta1, theta2, L1, L2, h);
    
    while norm(r,2)>tol
        J = solutions_Final_p2b_Jacobian(L2, theta2);
        sol = sol - J\r;
        x = sol(1); theta2 = sol(2);
        r = solutions_Final_p2b_Residual(x, theta1, theta2, L1, L2, h);
    end
    
end

function [r] = solutions_Final_p2b_Residual(x, theta1, theta2, L1, L2, h)
    r = zeros(2,1);
    r(1) = L1*cos(theta1) + L2*cos(theta2) - x;
    r(2) = L1*sin(theta1) + L2*sin(theta2) - h;
end

function [J] = solutions_Final_p2b_Jacobian(L2, theta2)
    
    J = zeros(2,2);
    J(1,1) = -1;
    J(1,2) = -L2*sin(theta2);
    J(2,2) = L2*cos(theta2);

end

