% FILE: h091.m
% uses myf1d.m
%
% myf(x,y) = myf(x,y,c) = c(1)*x^2 + c(2)*y^2 + c(3)*x*y + c(4)*x
% int myf(x,0) = (c(1)/3)*x^3+(c(4)/2)*x^2+C
% int myf(2,y) = (c(2)/3)*y^3+c(3)*y^2+(4*c(1)+2*c(4))*y

function h091
n=2; % dimension
mypath=[0 0 ; 2 0 ; 2 1]; % integration path

myf2 = inline('c(1)*x(:,1).*x(:,1)+c(2)*x(:,2).*x(:,2)+c(3)*x(:,1).*x(:,2)+c(4)*x(:,1)','x','c');
myf=inline('sin(x(:,1).*x(:,2)) + exp(x(:,1)-x(:,2))','x','y'); 

% test program
error = 0;
for i=1:500
   c = 10*rand(1,4)-5*rand(1,4);
   true = ( 8/3*c(1)+2*c(4) ) + ( c(2)/3+c(3)+4*c(1)+2*c(4) );
   estimate = mylineint(myf, mypath,c);
   error = error + (true-estimate)^2;
end
error

% line integral
function y = mylineint(funfcn ,intpath,c)
[d1,d2]=size(intpath);
s=0;
for j=1:(d1-1)
   a=intpath(j,:);   b=intpath(j+1,:);
   s=s+quad('myf1d',0,norm(a-b),[],[],a,b,funfcn,c);
end
y = s;
% FILE: h091.m
