If you have the Symbolic Toolbox, you might try replacing normcdf by an erf function and see what you get.
Otherwise you could try something like this:
x = linspace(-5,5);
a = pi;
f = @(x) 2.*(1/sqrt(2*pi)).*exp(-0.5*x.^2).*normcdf(a.*x,0,1);
F = zeros(size(x));
for j=1:length(x)
F(j) = integral(f,-Inf,x(j));
end
plot(x,F)
For the inverse cdf, you might try interpolating or using fzero.
