Problem with dsolve: Unable to find explicit solution

Hi everyone,

I am trying to solve the following ODE:

consider that [3] = a, [4] = b and [5] = c

I am using the following code to solve this:

answer = dsolve('Da=-k1*a + km1*b - k3*a*b','Db=k1*a - km1*b - k3*a*b','Dc=k3*a*b','a(0)=0.26','b(0)=0','c(0)=0')

However when I run this code I get the following:

Warning: Unable to find explicit solution. 
> In dsolve (line 201)
  In Paper2_Figure4 (line 13) 
 
answer =
 
[ empty sym ]

I don't understand what I'm doing wrong, any help would be greatly appreciated! Thanks :)

Answer 1

[Walter Roberson]

I found a way to do this, but I will need to wait until I get to my desktop to write through all of the steps as simplify() gave some unexpected results that are not understandable as-is.

The approach is to notice that the last sub-expression of the first two equations is the same as the expression for Dc. The first two are like Da=P+Dc , Db=-P+Dc in form. Select down to those two with the sub-expression replaced by Dc and dsolve. You will get solutions for a and c.

Take the expression for c and substitute t=0 to get c(0). You will get an expression in C1, C2, and some other terms and you know that c(0)=0 so solve() the expression for C2. The result will involve b(0) but we know that b(0)=0 so subs() that to get C2.

Now subs that C2 into the solution for a. subs t=0, solve() that equal 0.26 for C1, subs b(0)=0 and you get out C1

Now go back to the solutions for a and c and substitute in C1 and C2 and b(0)=0 to get the final a and c. Do not simplify() c!!! -- you will get weird results if you do!

Those in hand, go back to the second equation giving the expression for Db and substitute a and c and derivative of c to get Db.

At this point I simplify() and the dB term disappeared, and it was then that I noticed unexpected functions in the output, which I traced back to my having simplify() after I subs C1 into c. So at the moment I do not have the final production for b, but it should be down to a single equation and b(0)=0 will clearly be important. This resolution will have to wait until I get to my desktop.

[BobbyJoe]

Walter, I genuinely can't thank you enough! I greatly appreciate your time and effort to answer my question. Take all the time you need to get to your desktop :) In the mean time I will follow your explanation and see if I can crack it!

[BobbyJoe]

So far I am unsucessful in getting Db, I think I will just wait until you send the code so I can see what I'm doing wrong. Appreciate it!

[Walter Roberson]

syms a(t) b(t) c(t)
syms k1 km1 k3 real
syms C1 C2

Da = diff(a);
Db = diff(b);
Dc = diff(c);

eqn1 = Da == -k1*a + km1*b - k3*a*b;
eqn2 = Db == k1*a - km1*b - k3*a*b;
eqn3 = Dc == k3*a*b;
ic0_a = 0.26;
ic0_b = 0;
ic0_c = 0;

E1 = subs(eqn1, rhs(eqn3), lhs(eqn3));
E2 = subs(eqn2, rhs(eqn3), lhs(eqn3));
sol = dsolve(E1, E2);

a_partial = sol.a;
c_partial = sol.c;

c0 = subs(c_partial, t, 0);
C2sol = solve(c0 == ic0_c, C2);
C2solb = subs(C2sol, b(0), ic0_b);

a_partial2 = subs(a_partial, C2, C2solb);
a0 = subs(a_partial2, t, 0);

C1sol = solve(a0 == ic0_a, C1);
C1solb = subs(C1sol, b(0), ic0_b);

a_full1 = simplify(subs(subs(subs(a_partial, C2, C2solb), C1, C1solb), b(0), ic0_b));
c_full1 = subs(subs(subs(c_partial, C2, C2solb), C1, C1solb), b(0), ic0_b);   %DO NOT SIMPLIFY

%cross_check
if subs(subs(a_full1, t, 0), b(0), ic0_b) ~= ic0_a
    disp('cross-check on a(0) failed!')
else
    disp('cross-check on a(0) worked!');
end
if subs(subs(c_full1, t, 0), b(0), ic0_b) ~= ic0_c
    disp('cross-check on c(0) failed!')
else
    disp('cross-check on c(0) worked!');
end

b_partial = subs(subs(eqn2, a, a_full1), c, c_full1);
b_full = simplify(dsolve(b_partial, b(0) == ic0_b));

a_full = simplify(subs(a_full1, b, b_full));
c_full = simplify(subs(c_full1, b, b_full));

fprintf('a(t) == %s\n', char(a_full));
fprintf('b(t) == %s\n', char(b_full));
fprintf('c(t) == %s\n', char(c_full));

However! In practice it gets stuck when substituting in b_full into c_full1 . b_full is an expression involving multiple uses of the WittakerM function, and calculating the final c_full requires integrating it in two different places. a_full has a simple expression plus b(t) so no hard calculation is needed like for c_full

b_full is

(exp(-2*k1*t)*(2197*k3^3*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 2197*k3^3*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 16900*k1*k3^2*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 16900*k1*k3^2*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 8450*k3^2*km1*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 8450*k3^2*km1*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 12675*k1*k3^2*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 16900*k1*k3^2*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 2197*k3^2*abs(k3)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 2197*k3^2*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 4225*k3^2*km1*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 4225*k3^2*abs(k1 - km1)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 16900*k1*k3^2*exp(2*k1*t)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 12675*k1*k3^2*exp(2*k1*t)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 113750*k1^2*k3*exp(2*k1*t)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 195000*k1^2*k3*exp(2*k1*t)*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 195000*k1^2*k3*exp(2*k1*t)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 113750*k1^2*k3*exp(2*k1*t)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 97500*k1^2*exp(2*k1*t)*abs(k3)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 97500*k1^2*exp(2*k1*t)*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 16250*k3*km1^2*exp(2*k1*t)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 16250*k3*km1^2*exp(2*k1*t)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 4225*k3^2*km1*exp(2*k1*t)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 32500*km1^2*exp(2*k1*t)*abs(k3)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 32500*km1^2*exp(2*k1*t)*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 4225*k3*abs(k3)*abs(k1 - km1)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 130000*k1^2*k3*exp(2*k1*t)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 130000*k1^2*k3*exp(2*k1*t)*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 81250*k1^2*exp(2*k1*t)*abs(k3)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 65000*k1^2*exp(2*k1*t)*abs(k3)*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 65000*k1^2*exp(2*k1*t)*abs(k3)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 81250*k1^2*exp(2*k1*t)*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 16250*km1^2*exp(2*k1*t)*abs(k3)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 16250*km1^2*exp(2*k1*t)*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 16900*k1*k3*abs(k3)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 16900*k1*k3*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 8450*k3*km1*abs(k3)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 8450*k3*km1*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 4225*k3^2*exp(2*k1*t)*abs(k1 - km1)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 16900*k1*k3^2*exp(2*k1*t)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 16900*k1*k3^2*exp(2*k1*t)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 162500*k1^2*k3*exp(2*k1*t)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 162500*k1^2*k3*exp(2*k1*t)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 32500*k3*km1^2*exp(2*k1*t)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 32500*k3*km1^2*exp(2*k1*t)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 8450*k3^2*km1*exp(2*k1*t)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 8450*k3^2*km1*exp(2*k1*t)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 65000*k1^2*exp(2*k1*t)*abs(k3)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 65000*k1^2*exp(2*k1*t)*abs(k3)*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 12675*k1*k3*abs(k3)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 16900*k1*k3*abs(k3)*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 4225*k3*km1*abs(k3)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 16900*k1*k3*exp(2*k1*t)*abs(k3)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 16900*k1*k3*exp(2*k1*t)*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 65000*k1*k3*km1*exp(2*k1*t)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 65000*k1*k3*km1*exp(2*k1*t)*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 65000*k1*k3*km1*exp(2*k1*t)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 65000*k1*k3*km1*exp(2*k1*t)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 130000*k1*km1*exp(2*k1*t)*abs(k3)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 130000*k1*km1*exp(2*k1*t)*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 8450*k3*km1*exp(2*k1*t)*abs(k3)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 8450*k3*km1*exp(2*k1*t)*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 16900*k1*k3*exp(2*k1*t)*abs(k3)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 12675*k1*k3*exp(2*k1*t)*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 97500*k1*km1*exp(2*k1*t)*abs(k3)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 65000*k1*km1*exp(2*k1*t)*abs(k3)*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 65000*k1*km1*exp(2*k1*t)*abs(k3)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 97500*k1*km1*exp(2*k1*t)*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 4225*k3*km1*exp(2*k1*t)*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 48750*k1*k3*exp(2*k1*t)*abs(k1 - km1)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 48750*k1*k3*exp(2*k1*t)*abs(k1 - km1)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 16250*k3*km1*exp(2*k1*t)*abs(k1 - km1)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 16250*k3*km1*exp(2*k1*t)*abs(k1 - km1)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 130000*k1*k3*km1*exp(2*k1*t)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 130000*k1*k3*km1*exp(2*k1*t)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 65000*k1*km1*exp(2*k1*t)*abs(k3)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 65000*k1*km1*exp(2*k1*t)*abs(k3)*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 65000*k1*k3*exp(2*k1*t)*abs(k1 - km1)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 65000*k1*k3*exp(2*k1*t)*abs(k1 - km1)*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 16250*k1*exp(2*k1*t)*abs(k3)*abs(k1 - km1)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 16250*k1*exp(2*k1*t)*abs(k3)*abs(k1 - km1)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 4225*k3*exp(2*k1*t)*abs(k3)*abs(k1 - km1)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 16250*km1*exp(2*k1*t)*abs(k3)*abs(k1 - km1)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 16250*km1*exp(2*k1*t)*abs(k3)*abs(k1 - km1)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))))/(650*abs(k3)*(13*k3^2*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 13*k3^2*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 50*k1*k3*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 150*k1*abs(k3)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 150*k1*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 13*k3*abs(k3)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 13*k3*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 50*k3*km1*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 50*km1*abs(k3)*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 50*km1*abs(k3)*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 100*k1*abs(k3)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) + 200*k1*abs(k3)*whittakerW((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1)) - 50*k1*k3*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 50*k1*k3*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) + 50*k3*km1*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 50*k3*km1*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))*whittakerM(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1)) - 50*abs(k3)*abs(k1 - km1)*whittakerM((k3*km1 - k1*k3 + 4*k1*abs(k3))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*abs(k3))/(100*k1))*whittakerW(-(k3*(k1 - km1))/(4*k1*abs(k3)), -abs(k1 - km1)/(4*k1), -(13*exp(-2*k1*t)*abs(k3))/(100*k1))))

I do not know yet if there is a way to get MATLAB to rewrite that into simpler terms.

I recommend you put in a breakpoint on the a_full assignment so that you can observe that everything up to there has worked okay, that the dsolve calculation for b_full manages to work to form a closed-form expression. Just not necessarily one that integrates nicely in the context of finally resolving c.

[Walter Roberson]

By the way, Maple (a major symbolic product) does not have any more success in trying to solve the system without a lot of guideance.

[BobbyJoe]

Thank you so much Walter!! Really appreciate it!! I am fairly new to MATLAB so when I run this code, it comes with an error:

Error using solve>getEqns (line 404)
Input argument contains an empty equation or variable.

Error in solve (line 226)
[eqns,vars,options] = getEqns(varargin{:});

Error in test (line 34)
C1sol = solve(a0 == ic0_a, C1);
 

I'm not sure what I should do at this point :( But I really am grateful for all your help!

[Walter Roberson]

Ah, I assumed here that MATLAB would name the constants of integration C1 and C2, and although I can fetch the actual names using symvar, it would take me a bit more work to distinguish between them.

If I recall correctly, one of the two partials involves only the second constant (C2 here) and the other one involves both of them, so

C2var = setdiff(symvar(c_partial), [k1 km1 k3]);
C1var = setdiff(symvar(a_partial), [k1 km1 k3 C2var]);

However, I worry that your system might be returning empty from the dsolve() . Could you remind us which MATLAB release you are using?