% Example: Falling Ball with 2nd Order Runge Kutta
% Author: F. P. Schloerb
clear
g = 9.8;       % gravitational acceleration
X0 = 0;        % initial values
V0 = 0;
dt = 2;        % step size
t = 0:dt:100;  % array of times for calculation

X = zeros(length(t),1,'double');   % define solution arrays
V = zeros(length(t),1,'double');
XEuler = zeros(length(t),1,'double');

X(1) = X0;
XEuler(1) = X0;
V(1) = V0;
for i=1:length(t)-1
    % this is the 2o Runge-Kutta method
    Vmid = V(i) - g*dt/2;
    V(i+1) = V(i) - g*dt;
    X(i+1) = X(i) + Vmid*dt;

    % compute the Euler Method solution too for fun!
    XEuler(i+1) = XEuler(i) + V(i)*dt;
end

% make a nice graph
plot(t,XEuler,'md','MarkerSize',6)
hold on
plot(t,X,'ko','MarkerSize',6)
plot(t,X0+V0*t-0.5*g*t.*t,'k');
hold off
xlabel('t (s)');
ylabel('X (m)');
legend('Euler Method', '2nd Order RK', 'Exact')