xspan = [1 2];
alpha = 1;
beta = 2;
y10 = [alpha 0];
y20 = [0 1];

opts = odeset('RelTol',1e-8,'AbsTol',1e-8);
sol1 = ode45(@lineq1, xspan, y10);
sol2 = ode45(@lineq2, xspan, y20);

c2 = 1/70*(8 - 12* sin(log(2))- 4*cos(log(2)));
c1 = 11/10 - c2;

exact = @(x) c1*x + c2 ./ x.^2 - 3/10*sin(log(x)) - 1/10*cos(log(x));

x = linspace(1, 2, 50);
y1 = deval(sol1, x);
y2 = deval(sol2, x);
gamma = (beta - y1(1, end))/y2(1, end);
y = y1 + gamma*y2;

plot(x, y(1, :), 'b.');
hold on
plot(x, y1(1, :), '--');
plot(x, y2(1, :), '--');
plot(x, exact(x), 'r-')

title('Linear shooting method');
xlabel('x');
ylabel('y');
grid on
legend('Numerical solution', 'Aux 1', 'Aux 2', 'Analytical solution', 'Location', 'northwest')

function dydx = lineq1(x, y)
% LINEQ1 Evaluate the rhs of the Aux1 ODE
    dydx = [y(2); -2/x*y(2) + 2/x^2*y(1) + sin(log(x))/x^2];
end

function dydx = lineq2(x, y)
% LINEQ2 Evaluate the rhs of the Aux2 ODE
    dydx = [y(2); -2/x*y(2) + 2/x^2*y(1)];
end
