You were almost there. Hope this helps. The zero in the differential is a bit crude (just to keep the vectors the same length), but a small enough step size should give you are very accurate answer.
f = @(x) (2./sqrt(pi)).*integral(@(t) exp(-t.^2),0,x);
dx = 0.01;
x = (-5:dx:5)';
y = arrayfun(f,x);
dydx = [0;diff(y)./dx];
d2ydx2 = [0;diff(dydx)./dx];
[dy2d2x_min,minpos] = min(d2ydx2);
x_min = x(minpos)
figure;
plot(x,[y dydx d2ydx2],x_min,dy2d2x_min,'*')
legend('erf(x)','erf''(x)','erf''''(x)')
