Execute this code separately and copy the expressions obtained for az1 and bz1 into your code from above.
If you want more variability, also declare a0, r0, l and ta as syms. You'll get an expression for az1 and bz1 depending on t,x,y,z,a0,r0,l and ta. Make az1 and bz1 functions like
function value = AZ1(t,x,y,z,a0,r0,l,ta)
value = (expression obtained from the below code)
end
function value = BZ1(t,x,y,z,a0,r0,l,ta)
value = (expression obtained from the below code)
end
and call them in your code as
az1 = AZ1(t,x,y,z,a0,r0,l,ta)
bz1 = BZ1(t,x,y,z,a0,r0,l,ta)
syms t x y z
a0=10;
r0=70;
R=z* (1+ (r0^2/(2* z))^2);
zr=r0^2/2;
eps=r0/zr;
rho=sqrt((x^2+y^2)/r0^2);
f=sqrt(1+(2*z/r0^2)^2);
ph=(t-z)+2*atan(2*(z)/r0^2)-(x^2+y^2)/(2* R);
l=1;
ta = 200;
fr=(a0/f)*(sqrt(2)*rho/f)^l*exp(-rho^2/f^2)*exp(-((t-z)/ta)^2);
ax1=fr*cos(ph);
az1 = -diff(ax1,x)
bz1 = -diff(ax1,y)
%az1 = simplify(az1)
%bz1 = simplify(bz1)
az1 = matlabFunction(az1)
bz1 = matlabFunction(bz1)
