function [S,R]=simpson(f,a,b,n,M4)
% Syötteet:f   - integroitava funktio, tiedostossa vaikkapa intfun.m
%          a,b - v{lin p{{tepisteet
%            n - "pikkuv{lien lkm. parillinen, M4: f'''':n itseisarvon. max
%           Jos annat parittoman n:n, niin lasken (n+1):llä, ähäkutti! 
% Tulokset: S Simpsonarvo, R virhearvio

if nargin == 4
  M4=inf;         % Siitäs saat, jos olet niin typerä, että kutsut kahdella
end;              %  paluuargumentilla, vaikket anna M4:ää!
if nargin == 3    % Jos ei anneta n:ää, käytetään arvoa n=10.
  n=10;
end;
n=2*ceil(n/2);  % Parittomuuskorjaus
h=(b-a)/n;
w=2*ones(1,n-1);w(1:2:n-1)=4*ones(1,n/2); 
w=[1,w,1];
x=a:h:b;y=feval(f,x);
S=(h/3)*w*y';
if nargout==2
  R=(M4*(b-a)*h^4)/180;
end;
