clear
clc
close all

%% EDO
% y'' + y = g(t)
% y(0) = 0
% y'(0) = 1

%% Funcion g(t) con Heaviside

g = @(t)double (heaviside(t-pi) ...
              - heaviside(t-2*pi));

%% Funcion para ode45

fun = @(t,y) [y(2); -y(1) + g(t)];

%% Condiciones iniciales

y0 = [0 1];

%% Resolver

[t,y] = ode45(fun,[0 4*pi],y0);

%% Graficar

figure

plot(t,y(:,1),'LineWidth',2)

grid on
xlabel('t')
ylabel('y(t)')

title('y'''' + y = g(t)')

xline(pi,'--','t = \pi')
xline(2*pi,'--','t = 2\pi')

legend('Solucion ode45')