title: Is it possible to solve a system of PDEs using sympy in Python?
body: <p>Assume that you have a system of coupled PDEs, such as</p>
  
  <blockquote>
    <p>1st PDE in F(a,b)</p>
    
    <p>2nd PDE in F(a,b)</p>
  </blockquote>
  
  <p>Following code is able to solve each one of the PDEs separately:</p>
  
  <pre><code>import numpy as np
  import sympy as sp
  
  # definition of variables
  a, b = sp.symbols('a b')
  f = sp.Function('f')
  F = f(a, b)
  Fda = F.diff(a)
  Fdb = F.diff(b)
  
  # definition of PDEs
  eq1 = Fda - 2
  eq2 = Fda + Fdb + 2
  
  # solution of separated PDEs
  sp.pprint(sp.pdsolve(eq1))
  sp.pprint(sp.pdsolve(eq2))
  </code></pre>
  
  <p>Is it possible to solve the system of PDEs? The syntax would be something like <code>sp.pprint(sp.pdsolve([eq1, eq2]))</code>. I have tried <code>[eq1, eq2]</code>, <code>{eq1, eq2}</code>, <code>np.array([eq1, eq2])</code> etc. I had a look at <code>help(sp.pdsolve)</code> and <code>help(sp.pde)</code>, but have not yet found a solution.</p>