function [theta, omega, t] = solvependulum(L, tf, dt, theta0, omega0)
%SOLVEPENDULUM Numerical solution for simple pendulum motion
%   [THETA, OMEGA, T] = SOLVEPENDULUM(TF, DT, THETA0, OMEGA0)
%   solves for the angular position THETA and angular velocity
%   OMEGA of a simple pendulum having length L,
%   initial angular position THETA0 and initial angular
%   velocity OMEGA0 over the time period [0, TF] in time
%   steps of DT.


t = 0:dt:tf;
% make sure tf is at the end of t
if t(end) ~= tf
    t = [t tf];
end
n = length(t);

% estimated angular position and velocity
theta = zeros(size(t));
omega = zeros(size(t));
theta(1) = theta0;
omega(1) = omega0;
for i = 2:n
    theta(i) = theta(i - 1) + omega(i - 1) * dt;
    omega(i) = omega(i - 1) - 9.81 / L * sin(theta(i - 1)) * dt;
end
