Since x has t in the denominator. t=0 will cause problems, so start at t = 0.02 instead, perhaps! This would require a change in y to get the same number of elements as t.
You don't need x as a function. Replace with
x = y./((4*t*0.02/1000).^0.5);
then calculate your error function in a loop:
for i = 1:numel(t)
erf(i) = (2/pi^(0.5))*integral(fun,0,x(i));
end
Best not to use erf and erfc as they are functions already built-in to MATLAB. Better as, say, erF and erFc.
