
%  lab07 script to start with, integration
%

clear;

%  define the function to be integrated, and the bounds
%  anonymous function, note that x is a dummy variable
F_to_integrate = @(x) x*x*exp( -(x/3)^2 );  
%   integrate from a to b
a=0; b=10;  

%  fplot(F_to_integrate,[a,b]);  return; %.... plot the function 
%     (MATLAB complains, but it works)

%   another way to plot the function, no complaints, also logy
% xx=a:0.02:b; for i=1:length(xx); yy(i)=F_to_integrate(xx(i)); end; plot(xx,yy); return;

%  using your own function 'myint', integrate the function using different
%  methods and step sizes, and make a nice table, with each column
%  for the method and the rows for N
%  You should see how the different methods converge to the correct result

%  print a header line
fprintf('      N     Riemann-middle  Trapezoidal     Simpson-1/3 \n');
%  the output should all look like this:
   %   N     Riemann-middle  Trapezoidal     Simpson-1/3 
   %   5     4.001975562     0.116117360     4.447759508 
   %  25     1.730008408     1.696881107     1.728707786 
   %  45     1.713356227     1.713321689     1.717571040 
%  the numbers should line up nicely
%  use fprintf with proper formatting

for N=[5:20:900]   %  loop over the number of intervals
    %  now loop over the methods (a nested loop!)
    %  also, it loops over an array of strings
    fprintf('  %5d ',N);
    for meth = [ "Riemann-middle" "Trapezoidal" "Simpson-1/3" ] %  <-- this is a string array
        I=myint(F_to_integrate,a,b,N,meth);
        %  print the result here w/formatting
    end;  
    fprintf('\n');
end


function R=myint(f,a1,a2,n,m)   %  <--- dummy arguments can be different from actual arguments!
    if m == 'Riemann-middle'
        h = (a2-a1)/n;  %  size of the interval
        R=0;
        %  calculate integral
        return;
    end
    
    if m == 'Trapezoidal'  %  see lec14, page 21
        h = (a2-a1)/n;     %  size of the interval
        R=0;  
        %  calculate integral
        return;
    end
    
    if m == 'Simpson-1/3'  %  see lec14, page 36
        nn=rem(n,2);    %    the remainder
        n=n+nn;         %    make sure the number of points is odd, beware we start at i=0
        h = (a2-a1)/n;  %  size of the interval
        R=0;
        %  calculate integral
        return;
    end
    %  we shold never get here
    fprintf('error in function myint: no valid method given\n'); return;
end
