Plot Smile Volatility using Hull-White formula to Ornstein–Uhlenbeck process (Depend of Time and K)

I have the Hull-White formula:

if true
  % code
endif true
  % clear all
clc
theta = 1; a = 1.2; b = 0.3; % Vues in OU process
S0 = 100; r = 0; N =100; M = 1; T = 1; % Values in B-S formula
dt = T/N; 
sigma_path = zeros(M,N);
Xsquare = zeros(M,N);
tic
for j=1:M
  X = zeros(M,1); X(:,1) = 0;
  sigma = zeros(M,1); sigma(:,1) = 0;
  for k=1:N
    X(j,k+1) = X(j,k) + theta.*(a-X(j,k)).*dt + b.*sqrt(dt).*randn; %OU Process
    Xsquare(j,k) = X(j,k).*2; %square of OU Process
    sigma(j,k+1) = sigma(j,k) + Xsquare(j,k); 
    sigma_path(j,k) = sigma(j,k)./k; % Integral
    for i=1:80
      K(i) = 80 + (i./2);
      Time(k) = k./N;
      V(j,k,i) = BS_price(S0,K(i),r,Time(k),sqrt(sigma_path(j,k)));
      VV(k,i) = sum(V(j,k,i))./j; % average
      ImpVol(k,i) = BS_impv(S0,K(i),r, Time(k),VV(k,i));
    end
  end
end
toc
[mK,mT] = meshgrid(K,Time);
mesh(mK,mT,ImpVol)
end