function [T, U] = beam()
%BEAM Resolves the differential equation of a spring attached a beam
%   Returns the angle the beam makes with the vertical, and the angular
%   speed for 100 time-steps between 0 and 15.
%   Detailed explanation of the problem can be found in this document:
%   http://perso.uclouvain.be/vincent.legat/teaching/bac-q3/data/probleme-problem1415-5.pdf
%   (in french)

Tstart = 0;
Tend   = 15;

n = 1000; % Nombre de pas de temps
h = (Tend-Tstart)/n; % Taille d'un interval
T = linspace(Tstart,Tend,n+1); % Chaque pas de temps t
U = zeros(n+1,2); % Contient (theta theta') pour chaque temps t

% Condition initiale: quand t = 0
U(1, 1) = 0; % Theta = 0
U(1, 2) = 0; % Theta' = 0

% Application de Runge-Kutta d'ordre 4 sur chaque pas temps t
for i=1:n
    K1 = f(T(i),     U(i,:)        );
    K2 = f(T(i)+h/2, U(i,:)+h*K1/2 );
    K3 = f(T(i)+h/2, U(i,:)+h*K2/2 );
    K4 = f(T(i)+h,   U(i,:)+h*K3   );
    U(i+1,:) = U(i,:) + h*(K1+2*K2+2*K3+K4)/6;
end

end


function dudx = f(t,u)
dudx = u;

% Constante du problème
L = 1; % [m]
k = 140; % [N/m]
C = 10; % [Nms]
m = 10; % [kg]
g = 9.81; % [m/s^2]

sinTheta = sin(u(1));

dudx(1) = u(2);
dudx(2) = (L*(0.1-sinTheta*k*L) + m*g*L*sinTheta - C*u(2)) / (m*(2*L)^2/3);
end