% NymFys - Module 4.1 - Euler's Method - Part b
% Example 2 - Particle in free fall
% 17 October 2012

clear all;

format long % increases number of digits displayed

G=6.67300e-11; % gravitational constant
M=5.97219e24;  % mass of Earth
D=0.0025;      % drag coefficient   
m=1.0;         % mass of particle

N=100000; % number of steps
h=0.001; % step size

% initial values
t_0=0.0;
x_0=7.0e6;
v_0=0.0;

t=zeros(N+1,1);
x=zeros(N+1,1);
v=zeros(N+1,1);
t(1,1)=t_0;
x(1,1)=x_0;
v(1,1)=v_0;

t_old=t_0;
x_old=x_0;
v_old=v_0;

for n=1:N
    
    % Euler's method
    x_new=x_old+h*v_old; 
    v_new=v_old+h*((D/m)*v_old^2-G*m*M/x_old^2);
    
    t(n+1,1)=t_old+h;
    x(n+1,1)=x_new;
    v(n+1,1)=v_new;
    
    t_old=t_old+h;
    x_old=x_new;
    v_old=v_new;
    
end;

figure(1);
plot(t,x); % plotting the position vs. time: x(t)

figure(2);
plot(t,v); % plotting the velocity vs. time: v(t)

x_old % the value of x(t) at the last grid point t_N
v_old % the value of v(t) at the last grid point t_N
