function[Int]=intsp(N,xp)
%computes the convolution of f and g at points given by the matix xp on
%the interval [-1, 1] using N nodes gaussian quadriature

f=@(y) sin(y); %defines the funtion f
g=@(y) sin(exp(y)); %defines the function g
n=numel(xp); %finds the number of elements of xp
x=reshape(xp,n,1); %reshapes xp into a column vector
[z,w]=gauss(N); %finds the gaussian quadriture nodes z and weights w

%preallocate matrices for speed
t=zeros(n, N); 
Int=zeros(size(x)); 

for i=1:n
    t(i,:)=(f(x(i)-z)).*g(z);
    %computes values of the integrand at points given by xp and the quadrature nodes
       
end

Intp=t*w'; %sum_j(w_j f(xp-z_j)g(z_j) ie computes the integral using gaussian quad.
Int=reshape(Intp,size(xp)); %makes the integral the same shape as xp