%% Meteor closes point of approach
% Lewis 1.2-3

syms x1 x2 a b p1 p2 lambda

assume(x1>0)
assume(x2>0)

f = (x1-p1)^2 + (x2-p2)^2
h = x1^2/a^2 - x2^2/b^2 - 1

L = f + lambda*h

gL = gradient(L,[x1,x2,lambda])

sol = solve(gL,[x1,x2,lambda])

sol.x1
sol.x2
sol.lambda

%
% Even specifying the parameters a and b and position (p1,p2) will not help

gLs = subs(gL,[a,b,p1,p2],[1,2,3,4])

sol = solve(gLs,[x1,x2,lambda])

sol.x1
sol.x2
sol.lambda