tt=0;
close all;

M0=[1e16 3e14 1e14]';
X0=	[0 0;
	100 50;
	200 50]';

function retval = pot(x, t)

	M0=[1e16 3e14 1e14]';
	X0=	[0 0;
		100 50;
		200 50]';

	G=6.67e-11;

	retval=[0 0]';
	for l=1:columns(X0)
		retval+=G*M0(l)/norm(x-X0(:,l))^2*(X0(:,l)-x)/norm(x-X0(:,l));
	end
end

#x''=f(x,t), x'(t0)=v_0, x(t0)=x_0; здесь x - это радиус-вектор ракеты.
#
#раскладываем в систему уравнений:
#
#w: это [x, u, y, v]
#x'=u; u'=производная f по x;
#y'=v; v'=производная f по y;

function retval = f(w, t)
	retval=[];

	retval(2)=pot(w([1 3]), t)(1);
	retval(4)=pot(w([1 3]), t)(2);
	retval(1)=w(2);
	retval(3)=w(4);
end


#начальная скорость
v_0=[0 30]';

#начальное положение
x_0=[200 0]';

w0=[x_0(1) v_0(1) x_0(2) v_0(2)]';

t=0:0.005:30;

z=lsode("f", w0, t);

hold on;
plot(z(:,1), z(:,3));
plot(X0(1,:), X0(2,:), 'bo');
plot(x_0(1), x_0(2), 'r^');