Why isn't my animation of time varying isosurfaces animating, it seems to only capture one frame.

My aim is to create an animation of time varying isosurfaces corresponding to the poynting vector isosurfaces of a radiating electric dipole. Here is my code thus far

%1) Define grid

n = 100;
rmax = 2400;
prex = linspace(-rmax,rmax,n);
prey = linspace(-rmax,rmax,n);
prez = linspace(-rmax,rmax,n);

x = (nonzeros(prex))';
y = (nonzeros(prey))';
z = (nonzeros(prez))';

[X,Y,Z] = meshgrid(x,y,z);


%2) Define dipole moment magnitude, dipole moment and other constants

pmag = 0.001*(x(n-1)-x(n-2));
dir = [0,0,1];
p = pmag*dir;
 
k = 0.001*5*pi;
epsilon = 8.85*(10^(-12));
c = 3*((10)^(8));
c1 = (4*pi)^(-1);
c2= (4*pi*epsilon)^(-1);

%3) Probably should create an "r" vector/array and related r vector/array

r = sqrt(X.^2 + Y.^2 + Z.^2);
rnegative1 = r.^(-1);
rnegative2 = r.^(-2);
rnegative3 = r.^(-3);
rnegative4 = r.^(-4);
rnegative5 = r.^(-5);

%4) Need to define more constants

mu = (4*pi)*((10)^(-7));
omega = k*sqrt(mu*epsilon);

%5) Need to create a vector of times for which we want to create an animation, say one period

T = 2*pi*((omega)^(-1));
t = linspace(0,T,n);

%6) Define spatial dependence of H 

Hx = c*c1*exp(1i*k*r).*rnegative4.*((k^2)*(r.^2) + 1i*k*r).*(p(3)*Y - p(2)*Z);
Hy = c*c1*exp(1i*k*r).*rnegative4.*((k^2)*(r.^2) + 1i*k*r).*(p(1)*Z - p(3)*X);
Hz = c*c1*exp(1i*k*r).*rnegative4.*((k^2)*(r.^2) + 1i*k*r).*(p(2)*X - p(1)*Y);

%7) Define spatial dependence of E1

E1x = (k^2)*c2*exp(1i*k*r).*rnegative3.*(p(1)*p(3)*Z - (p(3)^2)*X - (p(2)^2)*X + p(1)*p(2)*Y);
E1y = (k^2)*c2*exp(1i*k*r).*rnegative3.*(p(1)*p(2)*X - (p(1)^2)*Y - (p(3)^2)*Y + p(2)*p(3)*Z);
E1z = (k^2)*c2*exp(1i*k*r).*rnegative3.*(p(2)*p(3)*Y - (p(2)^2)*Z - (p(1)^2)*Z + p(1)*p(3)*X);

%8) Define spatial dependence of E2

E2x = 3*c2*exp(1i*k*r).*rnegative5.*(ones(size(r)) - 1i*k*r).*(p(1)*X.^2 + p(2)*X.*Y + p(3)*X.*Z);
E2y = 3*c2*exp(1i*k*r).*rnegative5.*(ones(size(r)) - 1i*k*r).*(p(1)*X.*Y + p(2)*Y.^2 + p(3)*Y.*Z);
E2z = 3*c2*exp(1i*k*r).*rnegative5.*(ones(size(r)) - 1i*k*r).*(p(1)*X.*Z + p(2)*Y.*Z + p(3)*Z.^2);


%9) Define spatial dependence of E3

E3x = -c2*rnegative3.*exp(1i*k*r).*(ones(size(r)) - 1i*k*r).*rnegative3*p(1);
E3y = -c2*rnegative3.*exp(1i*k*r).*(ones(size(r)) - 1i*k*r).*rnegative3*p(2);
E3z = -c2*rnegative3.*exp(1i*k*r).*(ones(size(r)) - 1i*k*r).*rnegative3*p(3);


%10) Create the entire electric field

Ex = E1x + E2x + E3x;
Ey = E1y + E2y + E3y;
Ez = E1z + E2z + E3z;


%11) Create the color scheme for level sets

Scontx = real(Ey).*real(Hz) - real(Ez).*real(Hy);
Sconty = real(Ez).*real(Hx) - real(Ex).*real(Hz);
Scontz = real(Ex).*real(Hy) - real(Ey).*real(Hx);


Scontmag = sqrt(Scontx.^2 + Sconty.^2 + Scontz.^2);


Srangemax = max(Scontmag,[],'all');
Srangemin = min(Scontmag,[],'all');
Srange = linspace(Srangemin,Srangemax,100);


Srange(Srange == min(Srange)) = [];


%12) Create a for loop animation.

for i1 = 1:size(t)
    clf
    tt=t(i1)*ones(size(x));
    
    Ext = Ex.*exp(-1i*omega*tt);
    Eyt = Ey.*exp(-1i*omega*tt);
    Ezt = Ez.*exp(-1i*omega*tt);
    
    Hxt = Hx.*exp(-1i*omega*tt);
    Hyt = Hy.*exp(-1i*omega*tt);
    Hzt = Hz.*exp(-1i*omega*tt);
    
    Sx2 = real(Eyt).*real(Hzt) - real(Ezt).*real(Hyt);
    Sy2 = real(Ezt).*real(Hxt) - real(Hxt).*real(Hzt);
    Sz2 = real(Ext).*real(Hyt) - real(Eyt).*real(Hxt);
    
    
    Smag2 = sqrt(Sx2.^2 + Sy2.^2 + Sz2.^2);
    
    % Now the difficult (or rather the awkward) part, creating contour levels through arrayfun
    
    arrayfun(@(LEVEL) isosurface(X,Y,Z,Smag2,LEVEL), Srange);
    xlabel('x')
    ylabel('y')
    zlabel('z')
    
    % force matlab
    
    movieVector(i1) = getframe;
    
end

myWriter = VideoWriter('ThreeDdipolenoslices','MPEG-4');
myWriter.FrameRate = 20;

% Open the Video Writer object write the movie and close the file

open(myWriter);
writeVideo(myWriter,movieVector);
close(myWriter);

Smag2 is changing and is defined within the for loop (for t = time) and yet it doesn't seem to be changing at all in the plot results. I get one still frame. Furthermore, for some reason it shrinks my axis from (-2400,2400) to (-200,200).

It has to be something to do with Smag2, I must have inadvertantly made it static. At the very least I would expect the isosurface to change colors even if it's not moving.

I ultimately want to take slices of the varying isosurfaces and create something like the following

https://www.youtube.com/watch?v=UOVwjKi4B6Y

but instead of manually computing each cross section I'd like to just take slices.

Any help is greatly appreciated!