S = 100;
K = 110;
T = 1;
sigma = 0.1;
r = 0.02;
d = 0.0001;

d_plus = ((r + 1/2 * sigma^2) * T - log(K/S))/ (sigma * sqrt(T) );
d_minus = ((r - 1/2 * sigma^2) * T - log(K/S))/ (sigma * sqrt(T) );

delta_numerical = 1/d * (callBS(S+d, K,T,sigma,r) - callBS(S, K,T,sigma,r))
delta_theoretical = normcdf(d_plus)

vega_numerical = 1/d * (callBS(S, K,T,sigma+d,r) - callBS(S, K,T,sigma,r))
vega_theoretical1 = K*exp(-r*T)*normpdf(d_minus)*sqrt(T)
vega_theoretical2 = S*normpdf(d_plus)*sqrt(T)

rho_numerical = 1/d * (callBS(S, K,T,sigma,r + d) - callBS(S, K,T,sigma,r))
rho_theoretical = K*T*exp(-r*T)*normcdf(d_minus)

Gamma_numerical = 1/d^2*(callBS(S+d, K,T,sigma,r) -  2* callBS(S, K,T,sigma,r) + callBS(S-d, K,T,sigma,r) )
Gamma_theoretical = normpdf(d_plus)/(sigma*sqrt(T)*S)