y = linspace(0,2,100);
t = linspace(0,5,100);
[y1,t1] = meshgrid(y,t);
u_exact = 1-erf(y1./sqrt(4*t1));
u_gauss = 1-erf_gaussian(y1./sqrt(4*t1),10);
figure(1)
surf(y1,t1,u_exact)
xlabel(‘y’)
ylabel(‘t’)
figure(2)
surf(y1,t1,u_gauss)
xlabel(‘y’)
ylabel(‘t’)
function val = erf_gaussian(x,k)
[w,xi] = gauss_legendre_points(k);
[size_x,size_y] = size(x);
val = zeros(size(x));
for i=1:size_x
for j=1:size_y
f = @(t) exp(-(t.^2))/sqrt(pi);
t = 0.5*(1+xi)*x(i,j)-0.5*(1-xi)*x(i,j);
val(i,j) = dot(w,f(t))*x(i,j);
end
end