How to solve this differential equation inside derivation are there

function kk1
a0=10; r0=70; 
%main function to solve
[t,r]=ode45(@fn,[0 1200],[0.01 0.01 0.01 0.01 0.01 0.01])

plot(r(:,1),r(:,3))

    function dr = fn(t,r)
        % relation with actual variables
        x=r(1);y=r(2); z=r(3);px=r(4);py=r(5); pz=r(6);
        ta = 200;   
        zr=r0^2/2;
        rho=sqrt((x^2+y^2)/r0^2);
        eps=r0/zr;
        R=z* (1+ (r0^2/(2* z))^2); 
        f=sqrt(1+(2*z/r0^2)^2); 
        l=1;
        fr=(a0/f)*(sqrt(2)*rho/f)^l*exp(-rho^2/f^2)*exp(-((t-z)/ta)^2);
        ph=(t-z)+2*atan(2*(z)/r0^2)-(x^2+y^2)/(2* R));
        ax1=fr*cos(ph);
        ay1=0;
        az1=-d/dx(ax1);
        bx1=0;
        by1=ax1;
        bz1=- d/dy(ax1);
        g=sqrt(1+px^2+py^2+pz^2);
        dr=zeros(6,1);    
        dr(1)=px/g; dr(2)=py/g; dr(3)=pz/g; 
        dr(4)=-ax1*(1 - pz/g)-py/g* bz1; 
        
        dr(5)=-ay1+pz/g* ay1+px/g* bz1;
        dr(6)=-az1-px/g* ax1-py/g *ay1;
    end
end

here bz1=is the partial derivation w r to y of ax1 similar az1 is the partial diff w r to x with in the equation so how to solve this problem

Answer 1

[Torsten]

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)