With computer algebra systems like Mathematica, Maple, or sympy , it is all too tempting to use some single powerful command to solve a difficult problem. But without necessary care, you can quickly get into trouble! This is particularly true when trying to solve differential equations. Start a Jupyter notebook and consider the following simple nonlinear first-order equation using sympy :

We could solve it manually (and we will do so at the end of this article), but we could equally well just use the dsolve function in sympy :

This gives us the solution to the differential equation immediately. But really? We will see that this is not the case! And this is not sympy ’s fault. The same would happen in Mathematica or Maple. It’s the user’s fault because the user should always check the solution. Let’s see what’s wrong with this solution. Before we go deeper, let’s first fix the integration constant 𝐶1C1 by applying some initial condition, say 𝑦(0)=0:

So we can choose

to satisfy the initial condition.

We haven’t defined the symbol 𝐶1 in sympy . It was just given to us from the result of dsolve . So how can we substitute it by its actual value? We first need a handle on it. One way is to use free_symbols :

This gives us a set with all the free indices in an expression. Free indices are those which are not bound as summation indices, for instance. Sets are a bit impractical. You can iterate over it, but since there is no specified order, you cannot simply say: give me the first or only element. Instead, we can convert it to a list, which in our case would contain only one value, and then access it by index notation [0] :

Now we can substitute the constant:

Ok, so this is our solution from sympy . But is it right? Let’s take the derivative of this equation:

The left-hand side is the same as in our original differential equation, so substitute:

Let’s check if this equation is true for some specific values of 𝑥:

That’s ok.

Also ok.

Still ok.

Not ok! But it’s clear that our solution is wrong for some values of 𝑥. How come? Well, according to our differential equation, the slope of the solution must always be nonnegative, because the square root is defined as a positive quantity! Our “solution”, however, is a sine function, which obviously has a positive slope in some regions and a negative slope in other regions. So our solution cannot be right for all 𝑥! But where does this come from? The only way to find out is to solve the equation step by step instead of all in one step. So let’s do that now.

Start with our origin differential equation and put everything with 𝑦 on one side and everything with 𝑥 on the other side (this is called separation of variables). We get

Integrate the equation

and solve for y:

So this is the same solution that we got from dsolve . But where did we go wrong? We didn't, but we have made an implicit assumption. When we separated variables, putting everything with 𝑦 on one side and everything with 𝑥 on the other, we had to divide by

But if 𝑦²=1 for some 𝑥, we have a division by zero! Obviously, we have forgotten solutions, namely 𝑦=±1! In fact, the constant functions 𝑦=±1 satisfy the differential equation. But they don’t satisfy the initial condition 𝑦(0)=0. The only way to satisfy both the differential equation and the initial condition is to assemble the solution piecewise for different regions. One possible choice, and the only one which is continuous, is

Wrapping up: when solving differential equations with any computer algebra system, check the solution. It may be incomplete! One more reason why even nowadays it is imperative that you know your pen and paper methods.