clear all
clc
close all


deltat = 0.00001;
t = 0:deltat:4*pi;
N = length(t);

x = sin (t);

a1 = zeros(N, 1);
a2 = zeros(N, 1);
y  = zeros(N, 1);

y(1)  = 0;
a1(1) = 0;
a2(1) = 0;

n = 1;

while n <= N - 1

  a1(n + 1) = a2(n)*deltat + a1(n);
  a2(n + 1) = (-a2(n) - a1(n) + x(n) - y(n))*deltat + a2(n);
  y(n + 1)  = y(n) - deltat*y(n) + a1(n)*deltat;

  n = n +1
end

plot(a1); 
