clear, clc
format short g 
N = 100;
x = linspace(-pi,pi,N)'; % sine function over the interval [-π π]
x = (-pi:2*pi/(N-1):pi)';

y = sin(x); figure(1), plot(x,y) % plot sin(x) as a function of x
figure(1)
plot(x, y, 'b-', 'LineWidth', 2) % Blue solid line with increased width
grid on
xlabel('x (radians)', 'FontSize', 12)
ylabel('sin(x)', 'FontSize', 12)
title('Plot of sin(x) over [-\pi, \pi]', 'FontSize', 14)
ax = gca;
ax.FontSize = 12; % Increase font size for better readability
legend('sin(x)', 'Location', 'best')
%%
% Defining a function (II) 
sf = 0.1 ;
y = sin_scale(x,sf); % scale x by 0.1
figure(2), 
plot(x,y) 
%% 

sf=[0.1 0.5 2]; 	 	% Scale factor values in an array named sf 
N_sf = numel(sf); 	% Number of elements within sf (ML intrinsic % function) 

sin_array = zeros(N,N_sf); 
%Initialise an array of zeros to store data. % 

for ndx = 1: N_sf    
sin_array(:,ndx) = sin_scale(x,sf(ndx)); 
end 

plot(x, sin_array(:,1), 'r-', 'LineWidth', 2)   % Red solid line
hold on
plot(x, sin_array(:,2), 'b--', 'LineWidth', 2)  % Blue dashed line
plot(x, sin_array(:,3), 'g:', 'LineWidth', 2)   % Green dotted line
hold off

% Add Labels, Title, and Legend
xlabel('X values')
ylabel('Scaled Sine values')
title(['A Sine Function; Limits = [-\pi, \pi]; Scale Factors = ', num2str(sf)])
legend(arrayfun(@(s) sprintf('sf = %.1f', s), sf, 'UniformOutput', false), 'Location', 'Best')

% Improve Readability
grid on
set(gca, 'FontSize', 12)  % Increase font size for readability

%% Numerical Differentiation 
N = 100;

x = linspace(-pi,pi,N)';

function_y=cos(x);
% Plot function
figure(2)
plot(x,function_y)
xlabel('x')
ylabel('y')
title('y=cos(x)')
hold on
%% 
dy_by_dx = zeros(100, 1);

for i = 1:99
    delta_y = function_y(i+1) - function_y(i);
    delta_x = x(i+1) - x(i);
    
    % Store the result in dy_by_dx
    dy_by_dx(i) = delta_y / delta_x;
end

% For the last point, just copy the previous value to maintain the 100x1 size
dy_by_dx(100) = dy_by_dx(99);

plot(x,dy_by_dx,LineWidth=4,Color="b",Marker="_")
hold on
plot(x,-sin(x),Color="r",LineWidth=2)
hold off
% make this into a function 
%% 
f = @(x) x .* sin(x);
a= 0;
b = pi/2;
I = int_Sin(x);
plot(x, I)
TRAPEZ=trap1('int_Sin',a,b);
% per error
%% 
syms x
I = int(x * sin(x), x, 0, pi/2);
difference = ((TRAPEZ-1 )/1)* 100
disp(['The difference is for trapezoidal Method: ', num2str(difference), '%']);
%% Composite Trapexium Rule 
Composite = Comp(f, a, b, N); %function
% diff needs to be worked out 
%%
intervals = [5, 50, 100];%  intervals
[results_trap] = evaluate_function(intervals); %  function
% Display results
disp('Trapezoidal Rule Results:');
disp(results_trap);
intervals = [5, 50, 100]; %intervals
results_trap = evaluate_function(intervals);%  function
analytical_result = log(2); % analytical result
% Prepare the table
comparison_table = table(intervals', results_trap', ...
    repmat(analytical_result, size(intervals')), ...
    'VariableNames', {'Intervals', 'Trapezoidal_Result', 'Analytical_Result'});

% Display the comparison table
disp('Comparison of Numerical Integration Results with Analytical Result:');
disp(comparison_table);

%% Number of sub intervals needed to achieve accuracy to 6 decimal places for composite trapezium rule
f = @(x) 1 ./ x;
a = 1;
b = 2;
tolerance = 1e-6;
n1 = 1;
while true
   Composite = Comp(f, a, b, n1);
    if abs( Composite - analytical_result) <= tolerance
        break;
    end
    n1 = n1 + 1; 
end
disp(['Number of subintervals (n) used for trapezium rule: ', num2str(n1)]);
%% Simpsons rule
a = 1;         % Lower limit
b = 2;         % Upper limit
n_values = [5,50,100];       % intervals

%  array to store results
S = zeros(size(n_values));

% Loop over each value of n and calculate the integral
for i = 1:length(n_values)
    n = n_values(i);
    S(i) = Simpson(f, a, b, n); % Call your Simpson function
end

analytical_result = log(2); % Analytical result for comparison
% Create a table to display the results
results_table = table(n_values', S', repmat(analytical_result, size(S')), ...
    'VariableNames', {'n', 'Numerical_Result', 'Analytical_Result'});
% Display the table
disp("Numerical Result and analytical result for simpson")
disp(results_table);

%%
f = @(x) 1 ./ x;  % Define the function
a = 1;            % Lower limit
b = 2;            % Upper limit
analytical_result = log(2);  % Analytical result  ln(2)
tolerance = 1e-6; % Desired accuracy to six decimal places
n = 2;            % Start 

% Iteratively increase n 
while true
    I = compositeSimpson(f, a, b, n);
    if abs(I - analytical_result) <= tolerance
        break;
    end
    n = n + 1; 
end
% Display the result and number of intervals
disp(['Approximate integral: ', num2str(I)]);
disp(['Number of subintervals (n) used: ', num2str(n)]);

%% 
function [y]=my_funct(x)
    y=x*sin(x);
end