Solve a separable differential equation A first order separable ordinary differential equation is (arguably) the easiest type of differential equation to solve, and takes the form of

$$N(y)\frac{dy}{dx} = M(x) \\
y(x_0) = y_0$$ 

For two given functions \$N(y)\$ and \$M(x)\$ and an initial condition \$y(x_0) = y_0\$.

Your task is to take \$N(y)\$, \$M(x)\$, \$x_0\$ and \$y_0\$ and output the solution to that differential equation.

How to solve a separable differential equation

The differential equation

$$N(y)\frac{dy}{dx} = M(x)$$

can be solved through the following method:

$$
\begin{align}
N(y)\frac{dy}{dx} &amp; = M(x) \\
\int{N(y)\frac{dy}{dx}dx} &amp; = \int{M(x)dx} \\
\int{N(y)dy} &amp; = \int{M(x)dx}
\end{align}
$$

Once the integrals have been calculated, we can remove one of the constants of integration, leaving the equation

$$f(y) = g(x) + c$$

where

$$
f(y) = \int{N(y)dy} \\
g(x) = \int{M(x)dx}
$$

After this, \$c\$ can be calculated by applying the initial condition:

$$
c = f(y_0) - g(x_0)
$$

and a final solution for \$f(y)\$ can be given.

Challenge

You are to take in four inputs:


The function \$N(y)\$
The function \$M(x)\$
\$x_0\$
\$y_0\$


Both \$x_0\$ and \$y_0\$ can be any real number, and can be taken as input in many reasonable manner (i.e. a list, tuple, two integers etc.). The functions, in order to avoid having to focus on parsing mathematical input, will be limited to a type of Laurent polynomials - i.e polynomials in the form

$$
p(x) = \alpha \cdot x^n + \beta \cdot x^{n-1} + \dots + \gamma + \delta \cdot x^{-2} + \epsilon \cdot x^{-3} + \dots + \zeta \cdot x^{-m+1} + \eta \cdot x^{-m} \\
\alpha, \beta, \dots, \eta \in \mathbb{R}, \:\: n, m \in \mathbb{N}
$$

Notice that we'll never have \$x^{-1}\$ in the polynomial, and that because of this, the integrals will always be Laurent polynomials of the same form (with \$x^{-1}\$ included after integrating).

With inputs in this format, there are multiple ways you could represent them for input. As with the mapping, any reasonable representation is acceptable, and if you aren't sure, ask. However, here are some that are definitely allowed (for the example \$p(x) = 3x^2 + x + 6 - 2x^{-2} + 7x^{-5}\$):


Nested list: [[3, 2], [1, 1], [6, 0], [-2, -2], [7, -5]]
A pair of lists of coefficients: [3, 1, 6] and [0, -2, 0, 0, 7]
A string representation: "3x^2 + x + 6 - 2x^-2 + 7x^-5"
etc.


Given these four inputs, you are to output the solution to the differential equation they represent.

In order to avoid having to algebraically manipulate your solutions to get into the form \$y(x) = ...\$, you can output two Laurent polynomials, in the same form as the input; one representing \$y\$ and the other \$x\$.

This is code-golf so the shortest code in bytes wins!

Test cases

Both the MathJax and text-based inputs will be included for each example, where the text based will use the nested list input format above. In addition, I will walk through the first example for clarity.

In:
  N = [[1, -2]]
  M = [[6, 1]]
  1 -&gt; 1/25 = 0.04
Out:
  y = [[-1, -1]]
  x = [[3, 2], [28, 0]]


$$
\begin{align}
\frac{1}{y^2}\frac{dy}{dx} &amp; = 6x \\
y(1) &amp; = \frac{1}{25} = 0.04 \\
\\
\int{N(y)dy} &amp; = \int{M(x)dx} \\
\int{\frac{1}{y^2}dy} &amp; = \int{6xdx} \\
-\frac{1}{y} &amp; = 3x^2 + c \\
-25 &amp; = c + 3 \implies c = -28 \\
-\frac{1}{y} &amp; = 3x^2 - 28
\end{align}
$$

In:
  N = [[2, 1], [-4, 0]]
  M = [[3, 2], [4, 1], [-4, 0]]
  1 -&gt; 3
Out:
  y = [[1, 2], [-4, 1]]
  x = [[1, 3], [2, 2], [-4, 1], [-2, 0]]


$$
N(y) = 2y - 4, \:\: M(x) = 3x^2 + 4x - 4 \\
y(1) = 3 \\
\:\\
y^2 - 4y = x^3 + 2x^2 - 4x + c \\
c = -2 \\
y^2 - 4y = x^3 + 2x^2 - 4x - 2 \\
$$

In:
  N = [[1, 0]]
  M = [[3, 2], [2, 0]]
  0 -&gt; 0
Out:
  y = [[1, 1]]
  x = [[1, 3], [2, 1]]


$$
N(y) = 1, \:\: M(x) = 3x^2 + 2 \\
y(0) = 0 \\
\:\\
y = x^3 + 2x + c \\
c = 0 \\
y = x^3 + 2x \\
$$
Python 2, 123 122 bytes



def f(p,P):R=[[[a/-~b,b+1]for a,b in F]for F in P];a,b=map(lambda F,x:sum(a*x**b for a,b in F),R,p);R[1]+=[a-b,0],;print R


Try it online!

Straight-forward implementation. A function that takes arguments in the form (y, x), (N, M), and prints the resulting 2 polynomials to STDOUT.