#範例8-7：一階聯立微分方程式，求特解
#eq1：y1'(t) + 3y1 = 5
#eq1：y2'(t) - 3y2 = 5
#起始條件：y1(0)=1, y2(0)=2
#關鍵方法：dsolve()本身就可以解聯立微分方程式組
#語法：dsolve(聯立方程式list, [y1,y2]list)
#邊界條件：ics={y1.subs(t,0):1, y2.subs(t,0):2}
from sympy import *
t = symbols('t')
y1 = Function('y1')(t)
y2 = Function('y2')(t)
ode1 = Eq(y1.diff(t) + 3*y1, 5)
ode2 = Eq(y2.diff(t) - 3*y2, 5)
ans = dsolve([ode1, ode2], [y1, y2], ics={y1.subs(t,0):1, y2.subs(t,0):2})
print('y1=', ans[0].rhs)
print('y2=', ans[1].rhs)

#畫圖
f1 =ans[0].rhs
f2 =ans[1].rhs
#print(f)
p1 = plot(f1, f2, (t,-1,0.5), show=False)
#設定顏色
p1[0].line_color = 'r'
p1[1].line_color = 'g'
#設定圖例legend
p1[0].label = '$y1$'
p1[1].label = '$y2$'
#要打開 True
p1.legend = True
p1.show()