func.m

function z=func(y,x)
%invariabler är en vektor y och en skalär x, utvariabel en vektor
%av samma storlek som y
z=exp((x^2)*(y.^3));


dubbint.m

function int=dubbint(n)
%beräknar dubbelintegralen av 'func' över 0<=y<=fi=(x+1)exp(-x/2),
%a<=x<=b, med 2n delintervall på x-axeln (Simpsons formel)
a=0;b=1;
h=(b-a)/(2*n);x=a:h:b;
w=[0.5 ones(1,2*n-1) 0.5];
v=2:2:2*n;w(v)=2*ones(1,n);
fi=(x+1).*exp(-x/2);
vint=[];
for k=1:2*n+1
vint=[vint w(k)*quad('func',0,fi(k),[],[],x(k))];
end
int=(2*h/3)*sum(vint);


fcn.m

function out=fcn(t,x)
%t skall vara en vektor, x en skalär
fi=(x+1)*exp(-x/2);
out=fi*func(t*fi,x);


dubbint2.m

function int=dubbint2(a,b)
%ger dubbelintegralen av funk(y,x) över 0<=y<=fi=(x+1)*exp(-x/2),
%a<=x<=b med funktionen dblquad. fi ges i hjälpfunktionen fcn.
int=dblquad('fcn',0,1,a,b)