Can a solution from dsolve be used to feed equations for another instance of dsolve?

I found myself trying to solve the set of differential equations associated with uranium decay, but the system of equations for the decay series quickly becomes large enough for the computer to take far to long to solve. I had a thought, can I solve 3 or 4 steps with one dsolve, then take the last solution and feed it into another dsolve. Even better if I could do so after substituting alot of the variables with numbers, thus reducing the variable count considerably. But I cant seem to get anywhere, I assume im either doing something wrong, or its just not possible. This project is based off the bateman equation example found here https://www.mathworks.com/matlabcentral/fileexchange/69106-radioactive-decay-bateman-equation

This way works, but is slow, and as I add terms gets slower. the 12 term one below didnt solve in 3 hrs, up to term 8 took 2.5 hrs or so.

syms N1(t) N2(t) N3(t) N4(t) N5(t) N6(t) N7(t) N7(t) N8(t) N9(t) N10(t) N11(t) N12(t) l1 l2 l3 l4 l5 l6 l7 l8 l9 l10 l11 l12 r34 r35 r910 r911 N0; % l = lambda, r_xy = branching ratio 
eq_1 = diff(N1(t),t) == -l1*N1(t); % : decay of U-238 
eq_2 = diff(N2(t),t) == -l2*N2(t) + l1*N1(t); % formation of Th-234 from U-238 : decay of Th-234 
eq_3 = diff(N3(t),t) == -l3*N3(t) + l2*N2(t); % formation of Pa-234m from Th-234 : decay of Pa-234m 
eq_4 = diff(N4(t),t) == -l4*N4(t) + r34*l3*N3(t); % formation of Pa-234 from Pa-234m : decay of Pa-234 Branch 1 00.16% 
eq_5 = diff(N5(t),t) == -l5*N5(t) + l4*N4(t) + r35*l3*N3(t); % formation of U-234 from Pa-234 & Pa-234m : decay of U-234 Branch 2 99.84% 
eq_6 = diff(N6(t),t) == -l6*N6(t) + l5*N5(t); % formation of Th-230 from U-234 : decay of Th-230 
eq_7 = diff(N7(t),t) == -l7*N7(t) + l6*N6(t); % formation of Ra-226 from Th-230 : decay of Ra-226 
eq_8 = diff(N8(t),t) == -l8*N8(t) + l7*N7(t); % formation of Rn-222 from Ra-226 : decay of Rn-222 
eq_9 = diff(N9(t),t) == -l9*N9(t) + l8*N8(t); % formation of Po-218 from Rn-222 : decay of Po-218 
eq_10 = diff(N10(t),t) == -l10*N10(t) + r910*l9*N9(t); % formation of At-218 from Po-218 : decay of Po-218 Branch 1 00.16% 
eq_11 = diff(N11(t),t) == -l11*N11(t) + r911*l9*N9(t); % formation of Pb-214 from Po-218 : decay of Po-218 Branch 2 99.84% 
eq_12 = diff(N12(t),t) == -l12*N12(t) + l10*N10(t) + l11*N11(t);% formation of Bi-214 from At-218 & Pb-214 : decay of Bi-214

% solve the system of differential equations 
sol = dsolve ([eq_1 eq_2 eq_3 eq_4 eq_5 eq_6 eq_7 eq_8 eq_9 eq_10 eq_11 eq_12, N1(0)==N0 N2(0)==0 N3(0)==0 N4(0)==0 N5(0)==0 N6(0)==0 N7(0)==0 N8(0)==0 N9(0)==0 N10(0)==0 N11(0)==0 N12(0)==0]);

% solutions of the differential equations 
N1 = sol.N1;
N2 = sol.N2; 
N3 = sol.N3; 
N4 = sol.N4; 
N5 = sol.N5; 
N6 = sol.N6; 
N7 = sol.N7; 
N8 = sol.N8; 
N9 = sol.N9; 
N10 = sol.N10; 
N11 = sol.N11; 
N12 = sol.N12;

This is what i was trying to figure out, obviously i used a lower number of equations in the intital attempt, just 3 for each.

syms N1(t) N2(t) N3(t) N4(t) N5(t) N6(t) l1 l2 l3 l4 l5 l6 r34 r35 N0; % decay of U-238 to formation/decay of Pa-234m equations 
eq_1 = diff(N1(t),t) == -l1*N1(t); % : decay of U-238 
eq_2 = diff(N2(t),t) == -l2*N2(t) + l1*N1(t); % formation of Th-234 from U-238 : decay of Th-234 
eq_3 = diff(N3(t),t) == -l3*N3(t) + l2*N2(t); % formation of Pa-234m from Th-234 : decay of Pa-234m

sol = dsolve([eq_1 eq_2 eq_3, N1(0)==N0 N2(0)==0 N3(0)==0])
%Decay of Pa-234m into Th-230 equations %Where N3 is the solution N3 to the previous set of dsolve note N3s %presence in eqs 4 and 5. I also tried to put N3(t) in there 
eq_4 = diff(N4(t),t) == -l4*N4(t) + r34*l3*N3; % formation of Pa-234 from Pa-234m : decay of Pa-234 Branch 1 00.16% 
eq_5 = diff(N5(t),t) == -l5*N5(t) + l4*N4(t) + r35*l3*N3; % formation of U-234 from Pa-234 & Pa-234m : decay of U-234 Branch 2 99.84% 
eq_6 = diff(N6(t),t) == -l6*N6(t) + l5*N5(t); % formation of Th-230 from U-234 : decay of Th-230

sol2 = dsolve([eq_4 eq_5 eq_6, N4(0)==0 N5(0)==0 N6(0)==0])

% solutions of the differential equations 
N1 = sol.N1; 
N2 = sol.N2; 
N3 = sol.N3; 
N4 = sol.N4; 
N5 = sol.N5; 
N6 = sol.N6;

In the end, I may just give up on the full attempt for the decay series, and do one for only the longer lived prodgeny, assuming the shorter ones are just in secular equlibriem. But I really dont want this to be the solution to my issue. Any help is appreciated.

[Walter Roberson]

In the top section you include N5a(0) but there is no N5a ?

[John H.]

Hey, sorry Walter, thats a remnant of a previous attempt to split it up. When I copied the code Into the web I missed it.

I'll remove it.

Answer 1

[Walter Roberson]

syms N1(t) N2(t) N3(t) N4(t) N5(t) N6(t) l1 l2 l3 l4 l5 l6 r34 r35 N0;
% decay of U-238 to formation/decay of Pa-234m equations
eq_1 = diff(N1(t),t) == -l1*N1(t);                        %                                   : decay of  U-238
eq_2 = diff(N2(t),t) == -l2*N2(t) + l1*N1(t);             % formation of Th-234  from  U-238  : decay of Th-234
eq_3 = diff(N3(t),t) == -l3*N3(t) + l2*N2(t);             % formation of Pa-234m from Th-234  : decay of Pa-234m

sol = dsolve([eq_1 eq_2 eq_3, N1(0)==N0 N2(0)==0 N3(0)==0])

sol = struct with fields:
    N2: (N0*l1*exp(-l2*t))/(l1 - l2) - (N0*l1*exp(-l1*t))/(l1 - l2)
    N1: N0*exp(-l1*t)
    N3: (N0*l1*l2*exp(-l1*t))/((l1 - l2)*(l1 - l3)) - (N0*l1*l2*exp(-l2*t))/((l1 - l2)*(l2 - l3)) + (N0*l1*l2*exp(-l3*t))/((l1 - l3)*(l2 - l3))

%Decay of Pa-234m into Th-230 equations
    %Where N3 is the solution N3 to the previous set of dsolve note N3s
    %presence in eqs 4 and 5. I also tried to put N3(t) in there

N3 = sol.N3;

eq_4 = diff(N4(t),t) == -l4*N4(t) + r34*l3*N3;            % formation of Pa-234  from Pa-234m : decay of Pa-234  Branch 1  00.16%
eq_5 = diff(N5(t),t) == -l5*N5(t) + l4*N4(t) + r35*l3*N3; % formation of  U-234 from Pa-234 & Pa-234m  : decay of U-234  Branch 2  99.84%
eq_6 = diff(N6(t),t) == -l6*N6(t) + l5*N5(t);             % formation of Th-230 from   U-234  : decay of  Th-230

sol2 = dsolve([eq_4 eq_5 eq_6, N4(0)==0 N5(0)==0 N6(0)==0])

sol2 = struct with fields:
    N5: (exp(-l5*t)*(l5 - l6)*((N0*l1*l2*l3*l5*(l4*r34 + l4*r35 - l5*r35))/((l1 - l5)*(l2 - l5)*(l3 - l5)*(l4 - l5)*(l5 - l6)) - (N0*l1*l2*l3*l5*exp(l5*t - l1*t)*(l4*r34 + l4*r35 - l5*r35))/((l1 - l2)*(l1 - l3)*(l1 - l5)*(l4 - l5)*(l5 - l6)) + (N0*l1*…
    N4: -(exp(-l4*t)*(l4 - l5)*(l4 - l6)*((N0*l1*l2*l3*l4*l5*r34)/((l1 - l4)*(l4*l5 + l4*l6 - l5*l6 - l4^2)*(l2*l3 - l2*l4 - l3*l4 + l4^2)) + (N0*l1*l2*l3*l4*l5*r34*exp(l4*t - l1*t))/((l1 - l2)*(l1 - l3)*(l1 - l4)*(l4 - l5)*(l4 - l6)) - (N0*l1*l2*l3*l…
    N6: exp(-l6*t)*((N0*l1*l2*l3*l5*(l4*r34 + l4*r35 - l6*r35))/((l1 - l6)*(l2 - l6)*(l3 - l6)*(l4 - l6)*(l5 - l6)) - (N0*l1*l2*l3*l5*exp(l6*t - l1*t)*(l4*r34 + l4*r35 - l6*r35))/((l1 - l2)*(l1 - l3)*(l1 - l6)*(l4 - l6)*(l5 - l6)) + (N0*l1*l2*l3*l5*ex…

% solutions of the differential equations
N1 = sol.N1

N1 = 

N2 = sol.N2

N2 = 

N3 = sol.N3

N3 = 

N4 = sol2.N4

N4 = 

N5 = sol2.N5

N5 = 

N6 = sol2.N6

N6 = 

[John H.]

Thank you for your efforts Walter, I believe I attempted to define N3 in the way you defined it. But I will attempt your code this evening and post the results.