%% we are solving y''=-g  i.e  free fall motion
y0=[500,0];  % we start from the hight of 100 m  and our initial velocity is 0
timespan=[0,10]; % free fall for duration of  10 seconds

% independent variable x is time in our case

Npoints=10;
[time,y] = odeeuler(@free_fall_forces, timespan, y0, Npoints);

g=9.8;
yanalytical=y0(1) - g/2*time.^2;
vanalytical=y0(2) - g*time;

subplot(2,1,1);
plot(time, y(:,1), '-', time, yanalytical, 'r-');
legend('numerical','analytical');
xlabel('time');
ylabel('y-position');
title('position vs time');

subplot(2,1,2);
plot(time, y(:,2), '-', time, vanalytical, 'r-');
legend('numerical','analytical');
xlabel('time');
ylabel('y-velocity');
title('velocity vs time');
