# Solve 2 differential equations with 2 vairbales

I'm trying to solve a system of 2 differential equations with 2 variables. I know how to solve one differential equation. However, my things got a little messy with 2 equations. Here is what I did.

solAnalytic = DSolve[{
th1''[s] == fLinTh[th1[s], th2[s]],
th2''[s] == f2LinTh[th1[s], th2[s]],
th1[0] == 0.4,
th2[0] == 0.4,
th1'[0] == 0,
th2'[0] == 0}, {th1, th2}, s];


Where

fLinTh[theta1_, theta2_] = ((-m1 - m2) theta1)/m1 + (m2 theta2)/m1
f2LinTh[theta1_, theta2_] = (2 (l1 m1 + l1 m2) theta1)/(2 l1 m1 + l1 m2 + l2 m2) + (2 (l1 m1 + l1 m2) theta2)/(2 l1 m1 + l1 m2 + l2 m2)


At the end I would like to have 2 functions th1[theta1,theta2] and th2[theta1,theta2]

The result is a bit weird. Furthermore, I'm not sure how to replace the result of solAnalytic for 2 functions. I mean for a single function I usually the following for a system of one ode and one variable.

solution = DSolve[{
th''[s] == fLin[th[s]],
th[0] == 0.4, th'[0] == 0}, th, s];
theta[t_] = th[t] /. solution[[1]];


Where

fLin[th[s]] = -th[s]


And

theta[6] = 0.096017


Any help will be appreciated. I'm a bit lost.

Thank you

Edit: After some trial and error, I found that even if the function is not v1 and v2 dependent, I need to add those argument, because the function was created with those argument earlier. Thus,

solAnalytic = DSolve[{
th1''[s] == fLinTh[th1[s], th2[s], v1, v2],
th2''[s] == f2LinTh[th1[s], th2[s], v1, v2],
th1[0] == 0.4,
th2[0] == 0.4,
th1'[0] == 0,
th2'[0] == 0}, {th1, th2}, s];


Where

fLinTh[th1[s], th2[s], v1, v2] = -((7 th1[s])/5) + (2 th2[s])/5
f2LinTh[th1[s], th2[s], v1, v2] = (14 th1[s])/15 + (14 th2[s])/15


However, the solution seems weird...

Sol[s_] = th1[s] /. solAnalytice[[1]];
Sol[t] = 4.5 RootSum[-47250 + 105 #1^2 + #1^4 &,
E^((t #1)/15)/(105 + 2 #1^2) &] +
0.05 RootSum[-47250 + 105 #1^2 + #1^4 &, (-210 E^((t #1)/15)
+E^((t #1)/15) #1^2)/(105 + 2 #1^2) &]


I'm not sure how to replace the result of solAnalytic for 2 functions

One way could be

ClearAll["Global*"]
fLinTh[theta1_, theta2_] = ((-m1 - m2) theta1)/m1 + (m2 theta2)/m1
f2LinTh[theta1_, theta2_] = (2 (l1 m1 + l1 m2) theta1)/(2 l1 m1 + l1 m2 + l2 m2) + (2 (l1 m1 + l1 m2) theta2)/(2 l1 m1 + l1 m2 + l2 m2)

ode1 = th1''[s] == fLinTh[th1[s], th2[s]]
ode2 = th2''[s] == f2LinTh[th1[s], th2[s]]
ic = {th1[0] == 4/10, th2[0] == 4/10, th1'[0] == 0, th2'[0] == 0}


And now

solAnalytic = DSolveValue[{ode1, ode2, ic}, {th1[s], th2[s]}, s];
th1Solution[s_, m1_, m2_, l1_, l2_] = First@solAnalytic;
th2Solution[s_, m1_, m2_, l1_, l2_] = Last@solAnalytic;


Then you can do

th1Solution[s, 1, 2, 3, 4]


th2Solution[s, 1, 2, 3, 4]


You can use Manipulate to change m1,m1,l1,l2 and plot the solutions also for some range of s. ps. better to use exact numbers with DSolve

• By using the method I wrote, I get (3. (-0.2 + 0.0333333 Cos[2 theta1 - 2 theta2] + s^2 Sin[theta1] + 0.166667 s^2 Sin[theta1 - 2 theta2]))/(-6. + Cos[2 theta1 - 2 theta2])  for theta1Double with m1 and m2 fixed. Which is more what I expect. However, for some reason theta1 and theta2 doesn't get the initiales conditions. Commented Dec 4, 2022 at 23:49
• However, for some reason theta1 and theta2 doesn't get the initiales conditions. I do not understand what you mean. You do have ic = {th1[0] == 4/10, th2[0] == 4/10, th1'[0] == 0, th2'[0] == 0} set there and the ic is used in  DSolveValue[{ode1, ode2, ic}, {th1[s], th2[s]}, s];` so it all look OK to me. Are you saying DSolve is not using the initial conditions in the solution? I am not sure I understand what you mean. Commented Dec 5, 2022 at 0:04
• I will edit my post to give an example for one ode and one variable. For this system the code work as expected. Commented Dec 5, 2022 at 0:23