function IV = IV_BS(Data, Spot, r, T, flag)
IV = zeros(size(Data(:,1)));
options = optimoptions('fsolve','FunctionTolerance',1e-10,'StepTolerance',1e-10,'OptimalityTolerance',1e-10,'MaxIterations',1000,'MaxFunctionEvaluations',1000,'Display','off');

% Computation of the implied volatilities
for i = 1:length(IV)
    sigma = fsolve(@(vol) objfun_BS(vol,Spot,Data(i,2),r,T,Data(i,1),flag),0.45,options);
    IV(i) = sigma;
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dp = objfun_BS(vol, Spot, Strike, r, T, Price, flag)
% Computation of the B&S Price
d1 = (log(Spot./Strike)+(r+vol.*vol/2).*T)./(vol*sqrt(T));
d2 = d1-vol.*sqrt(T);

if flag == 1      % 1 for Call
    nd1 = normcdf(d1,0,1); 
    nd2 = normcdf(d2,0,1);
    bs = Spot.*nd1-Strike.*exp(-r.*T).*nd2;
elseif flag == -1     % -1 for Put
    nd1 = normcdf(-d1,0,1); 
    nd2 = normcdf(-d2,0,1);
    bs = Strike.*exp(-r.*T).*nd2-Spot.*nd1;
end

% Difference between the B&S price and the actual price
dp = bs-Price;

end