how to solve this differential equation within acceptable time complexity , ode45or ode23 seems uncapable of solving this

clc, clear
u_init = 10;
v_init = 0;
r_init = 30;
y0 = [u_init; v_init; r_init];
sigma = 20;  
np = 10;     
tspan = [0 , 100]; 
[t, y] = ode45(@(t, y) myODEs(t, y, sigma, np), tspan, y0);
figure;
plot(t, y(:, 1), 'r', t, y(:, 2), 'g', t, y(:, 3), 'b', 'LineWidth', 2);
legend('u', 'v', 'r');
xlabel('Time');
ylabel('Values');
title('Numerical Solution using ode15s with Smaller Time Steps');



function dydt = myODEs(t, y, sigma, np)
% Extract the variables from the input vector y
u = y(1);
v = y(2);
r = y(3);
load motion_equation
% Define your system of differential equations
du_dt = (1099511627776*(24500*u^2 + 24500*v^2)*((539*r^2)/(1000*(u^2 + v^2)) - v^2/(25*(u^2 + v^2)) + (89*v^4)/(125*(u^2 + v^2)^2) + (7*r*v)/(500*(u^2 + v^2)) - 41/2000))/3595427212083331 - (489383537374431765*np^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(11664*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(2160*np) - 2931/10000))/942519671084372721664 + (3595648213920514*r*v)/3595427212083331 + (898850755706880*r^2)/3595427212083331 - (84250078478336*sin(sigma)*sin(sigma - (79*(u^2 + v^2)^(1/2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/(218*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115)^(1/2)))*((32512578371218140519*(u^2 + v^2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180)^2)/2814749767106560000 + (61894238684256165279*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115))/703687441776640000))/449428401510416375;
dv_dt = (3024811283085532177211850752*(24500*u^2 + 24500*v^2)*((581*r)/(1000*(u^2 + v^2)^(1/2)) - (59*v)/(200*(u^2 + v^2)^(1/2)) + (343*r^3)/(125*(u^2 + v^2)^(3/2)) - (841*v^3)/(500*(u^2 + v^2)^(3/2)) + (2653*r*v^2)/(1000*(u^2 + v^2)^(3/2)) - (18767*r^2*v)/(1000*(u^2 + v^2)^(3/2))))/9689824258347110843885758248353 - (19378432542558421171785885212899*r*u)/19379648516694221687771516496706 + (247074214383739837580574720*(171500*u^2 + 171500*v^2)*((343*r)/(1000*(u^2 + v^2)^(1/2)) + (127*v)/(1000*(u^2 + v^2)^(1/2)) + (1029*r^3)/(250*(u^2 + v^2)^(3/2)) + (3*v^3)/(100*(u^2 + v^2)^(3/2)) + (1029*r*v^2)/(500*(u^2 + v^2)^(3/2)) - (2401*r^2*v)/(1000*(u^2 + v^2)^(3/2))))/9689824258347110843885758248353 - (496069050426027277062743523328*cos(sigma)*sin(sigma - (79*(u^2 + v^2)^(1/2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/(218*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115)^(1/2)))*((32512578371218140519*(u^2 + v^2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180)^2)/2814749767106560000 + (61894238684256165279*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115))/703687441776640000))/1211228032293388855485719781044125 - (247074214383739837580574720*cos(sigma)*sin(sigma - (79*(u^2 + v^2)^(1/2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/(218*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115)^(1/2)))*((32512578371218140519*(u^2 + v^2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180)^2)/2814749767106560000 + (61894238684256165279*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115))/703687441776640000)*((613*sin(sigma)*sin(sigma - (79*(u^2 + v^2)^(1/2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/(218*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115)^(1/2)))*((32512578371218140519*(u^2 + v^2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180)^2)/2814749767106560000 + (61894238684256165279*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115))/703687441776640000))/1000 + 15834/15625))/9689824258347110843885758248353;
dr_dt = (8104034231786666672642850816*cos(sigma)*sin(sigma - (79*(u^2 + v^2)^(1/2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/(218*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115)^(1/2)))*((32512578371218140519*(u^2 + v^2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180)^2)/2814749767106560000 + (61894238684256165279*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115))/703687441776640000))/242245606458677771097143956208825 - (247074214383739837580574720*(24500*u^2 + 24500*v^2)*((581*r)/(1000*(u^2 + v^2)^(1/2)) - (59*v)/(200*(u^2 + v^2)^(1/2)) + (343*r^3)/(125*(u^2 + v^2)^(3/2)) - (841*v^3)/(500*(u^2 + v^2)^(3/2)) + (2653*r*v^2)/(1000*(u^2 + v^2)^(3/2)) - (18767*r^2*v)/(1000*(u^2 + v^2)^(3/2))))/9689824258347110843885758248353 - (988364255149402852704649216*(171500*u^2 + 171500*v^2)*((343*r)/(1000*(u^2 + v^2)^(1/2)) + (127*v)/(1000*(u^2 + v^2)^(1/2)) + (1029*r^3)/(250*(u^2 + v^2)^(3/2)) + (3*v^3)/(100*(u^2 + v^2)^(3/2)) + (1029*r*v^2)/(500*(u^2 + v^2)^(3/2)) - (2401*r^2*v)/(1000*(u^2 + v^2)^(3/2))))/9689824258347110843885758248353 - (49661917091137100458229760*r*u)/9689824258347110843885758248353 + (988364255149402852704649216*cos(sigma)*sin(sigma - (79*(u^2 + v^2)^(1/2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/(218*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115)^(1/2)))*((32512578371218140519*(u^2 + v^2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180)^2)/2814749767106560000 + (61894238684256165279*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115))/703687441776640000)*((613*sin(sigma)*sin(sigma - (79*(u^2 + v^2)^(1/2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/(218*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115)^(1/2)))*((32512578371218140519*(u^2 + v^2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180)^2)/2814749767106560000 + (61894238684256165279*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115))/703687441776640000))/1000 + 15834/15625))/9689824258347110843885758248353;
% Pack the derivatives into the output vector dydt
dydt = [du_dt; dv_dt; dr_dt];
end

[Star Strider]

The ‘motion_equation’ file seems to be missing.

That aside, it should be loaded before the ‘myODEs’ function, and passed as an argument rather than loading it inside the function each time the function is called.

[子青]

Sir, I've seen you many times, truly appreciate your selfless contribution. Just wanted you to know how grateful I am

[Star Strider]

Thank you!

Answer 1

[Sam Chak]

Hi @子青

I used ode15s, and it appears that the system becomes very stiff around the time 0.0015. The equations are quite long. Please check if the terms and operators are entered correctly.

clc, clear
u_init  = 10;
v_init  = 0;
r_init  = 30;
y0      = [u_init; v_init; r_init];
sigma   = 20;  
np      = 10;     
% tspan   = [0, 100]; 
tspan   = [0, 1]; 
[t, y]  = ode15s(@(t, y) myODEs(t, y, sigma, np), tspan, y0);

Warning: Failure at t=1.449308e-03. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.469447e-18) at time t.

figure;
plot(t, y(:, 1), 'r', t, y(:, 2), 'g', t, y(:, 3), 'b', 'LineWidth', 2);

Warning: Imaginary parts of complex X and/or Y arguments ignored.

legend('u', 'v', 'r');
xlabel('Time');
ylabel('Values');
title('Numerical Solution using ode15s with Smaller Time Steps');

function dydt = myODEs(t, y, sigma, np)
    % Extract the variables from the input vector y
    u = y(1);
    v = y(2);
    r = y(3);
    % load motion_equation
    
    % Define your system of differential equations
    du_dt = (1099511627776*(24500*u^2 + 24500*v^2)*((539*r^2)/(1000*(u^2 + v^2)) - v^2/(25*(u^2 + v^2)) + (89*v^4)/(125*(u^2 + v^2)^2) + (7*r*v)/(500*(u^2 + v^2)) - 41/2000))/3595427212083331 - (489383537374431765*np^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(11664*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(2160*np) - 2931/10000))/942519671084372721664 + (3595648213920514*r*v)/3595427212083331 + (898850755706880*r^2)/3595427212083331 - (84250078478336*sin(sigma)*sin(sigma - (79*(u^2 + v^2)^(1/2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/(218*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115)^(1/2)))*((32512578371218140519*(u^2 + v^2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180)^2)/2814749767106560000 + (61894238684256165279*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115))/703687441776640000))/449428401510416375;
    dv_dt = (3024811283085532177211850752*(24500*u^2 + 24500*v^2)*((581*r)/(1000*(u^2 + v^2)^(1/2)) - (59*v)/(200*(u^2 + v^2)^(1/2)) + (343*r^3)/(125*(u^2 + v^2)^(3/2)) - (841*v^3)/(500*(u^2 + v^2)^(3/2)) + (2653*r*v^2)/(1000*(u^2 + v^2)^(3/2)) - (18767*r^2*v)/(1000*(u^2 + v^2)^(3/2))))/9689824258347110843885758248353 - (19378432542558421171785885212899*r*u)/19379648516694221687771516496706 + (247074214383739837580574720*(171500*u^2 + 171500*v^2)*((343*r)/(1000*(u^2 + v^2)^(1/2)) + (127*v)/(1000*(u^2 + v^2)^(1/2)) + (1029*r^3)/(250*(u^2 + v^2)^(3/2)) + (3*v^3)/(100*(u^2 + v^2)^(3/2)) + (1029*r*v^2)/(500*(u^2 + v^2)^(3/2)) - (2401*r^2*v)/(1000*(u^2 + v^2)^(3/2))))/9689824258347110843885758248353 - (496069050426027277062743523328*cos(sigma)*sin(sigma - (79*(u^2 + v^2)^(1/2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/(218*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115)^(1/2)))*((32512578371218140519*(u^2 + v^2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180)^2)/2814749767106560000 + (61894238684256165279*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115))/703687441776640000))/1211228032293388855485719781044125 - (247074214383739837580574720*cos(sigma)*sin(sigma - (79*(u^2 + v^2)^(1/2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/(218*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115)^(1/2)))*((32512578371218140519*(u^2 + v^2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180)^2)/2814749767106560000 + (61894238684256165279*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115))/703687441776640000)*((613*sin(sigma)*sin(sigma - (79*(u^2 + v^2)^(1/2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/(218*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115)^(1/2)))*((32512578371218140519*(u^2 + v^2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180)^2)/2814749767106560000 + (61894238684256165279*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115))/703687441776640000))/1000 + 15834/15625))/9689824258347110843885758248353;
    dr_dt = (8104034231786666672642850816*cos(sigma)*sin(sigma - (79*(u^2 + v^2)^(1/2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/(218*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115)^(1/2)))*((32512578371218140519*(u^2 + v^2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180)^2)/2814749767106560000 + (61894238684256165279*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115))/703687441776640000))/242245606458677771097143956208825 - (247074214383739837580574720*(24500*u^2 + 24500*v^2)*((581*r)/(1000*(u^2 + v^2)^(1/2)) - (59*v)/(200*(u^2 + v^2)^(1/2)) + (343*r^3)/(125*(u^2 + v^2)^(3/2)) - (841*v^3)/(500*(u^2 + v^2)^(3/2)) + (2653*r*v^2)/(1000*(u^2 + v^2)^(3/2)) - (18767*r^2*v)/(1000*(u^2 + v^2)^(3/2))))/9689824258347110843885758248353 - (988364255149402852704649216*(171500*u^2 + 171500*v^2)*((343*r)/(1000*(u^2 + v^2)^(1/2)) + (127*v)/(1000*(u^2 + v^2)^(1/2)) + (1029*r^3)/(250*(u^2 + v^2)^(3/2)) + (3*v^3)/(100*(u^2 + v^2)^(3/2)) + (1029*r*v^2)/(500*(u^2 + v^2)^(3/2)) - (2401*r^2*v)/(1000*(u^2 + v^2)^(3/2))))/9689824258347110843885758248353 - (49661917091137100458229760*r*u)/9689824258347110843885758248353 + (988364255149402852704649216*cos(sigma)*sin(sigma - (79*(u^2 + v^2)^(1/2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/(218*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115)^(1/2)))*((32512578371218140519*(u^2 + v^2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180)^2)/2814749767106560000 + (61894238684256165279*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115))/703687441776640000)*((613*sin(sigma)*sin(sigma - (79*(u^2 + v^2)^(1/2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/(218*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115)^(1/2)))*((32512578371218140519*(u^2 + v^2)*((497*r)/(100*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180)^2)/2814749767106560000 + (61894238684256165279*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((72*((1 - (15625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2*((34625*u^2*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70)^2)/(1458*np^2) + (2753*u*(exp(-2*abs((7*r)/(2*(u^2 + v^2)^(1/2)) + (pi*atan(v/u))/180))/70 - 11/70))/(270*np) - 2931/1250))/(729*np^2*pi))^(1/2)/2 + 1/2)^2)/115 + 43/115))/703687441776640000))/1000 + 15834/15625))/9689824258347110843885758248353;
    
    % Pack the derivatives into the output vector dydt
    dydt = [du_dt; 
            dv_dt; 
            dr_dt];
end