% Octave code for the solution of Exercise 1
%
% Program starts
%
% clear the variables and close the windows
%
clear all
close all
%
% Exercise 1
%
% With Octave
printf('Hello, this is the beginning ! \n');
%
% With MatLab
% display('Hello, this is the beginning !');
%
%  Definition of the integration range (sampled)
%
dt = 0.2;
t = 0:dt:8;
%
% Acceleration
%
a = - t .^3 + 10 * t.^2 - 20 * t ;
%
% velocity (analytic solution)
% initial condition: if v = 0 then t = 0
%
v = - (1/4) * t.^4 + (10/3)*t.^3 - 10*t.^2;
%
% space (analytic solution)
%
s = - (1/20)* t.^5 + (5/6)*t.^4 - (10/3)*t.^3;
%
% constant null function
%
zero = s * 0;
%
% I open a window for a figure, now
figure(1)
plot(t,a,'-b','LineWidth',1);
ylabel(' Space = green, Velocity = red, Acceleration = blue');
xlabel(' Time in s ');  
title([' Analytic method (hey folks I can also write Greek letters \mu ) ']);
set(1,'Color','w');
hold on
plot(t,v,'-r','LineWidth',1);
hold on
plot(t,s,'-g','LineWidth',1);
hold on
plot(t,zero,'-.k','LineWidth',1);
hold off
%
%
%    BEZOUT FORMNULA
%
vB = zero;
termine = length(zero)-1;
for j = 1:1:termine 
    vB(j+1) = vB(j) + dt * ( (a(j+1) + a(j))/2 );
end;
%
%
%   CAVALIERI SIMPSON
%
vCS = zero;
vCS(2) = v(2);
for j = 1:2:termine 
    vCS(j+2) = vCS(j) + dt * ( a(j) + 4 * a(j+1) + a(j+2) ) / 3 ;
end;
for j = 2:2:(termine-1) 
    vCS(j+2) = vCS(j) + dt * ( a(j) + 4 * a(j+1) + a(j+2) ) / 3 ;  
end
%
figure(2)
plot(t,v,'-b','LineWidth',1);
ylabel(' Velocity (analytic) = blue, Bezout = red , Cavalieri Simpson = green');
xlabel(' Time in s ');  
title([' Bezout and Cavalieri Simpson (example of subscript B_3 ) ']);
set(2,'Color','w');
hold on
plot(t,vB,'-.r','LineWidth',1);
hold on
plot(t,vCS,'-.g','LineWidth',1);
hold on
plot(t,zero,'-.k','LineWidth',1);
hold off
%
%  Velocity errors (they depend on the dt accuracy !!!
%
erroriB = vB - v;
erroriCS = vCS - v;
figure(3)
ylabel(' Bezout errors = red , Cavalieri Simpson errors = green');
xlabel(' Time in s ');  
title([' Bezout and Cavalieri Simpson Errors (for dt = ',num2str(dt), ') ']);
set(3,'Color','w');
hold on
plot(t,erroriB,'-.r','LineWidth',1);
hold on
plot(t,erroriCS,'-.g','LineWidth',1);
hold off
%
% Space integration
% Bezout and Cavalieri Simpson
%
%  .... Now it is your turn ! ............
%
%
% with Matlab
% display(' I have done ! ')
% display(' Good bye class ! ')
% with Octave
printf('I have done! \n');
printf('Good bye class ! \n');