Я не физик, поэтому не могу оценить вопрос ТС с этой точки зрения. Наверно в этом случае, он бы задал вопрос в другом разделе. Но вопрос был задан в математическом разделе и автор темы надеялся получить именно такую консультацию. Его интересовали численные методы решения дифференциального уравнения первого порядка с неявно выраженной производной для анализа фазовых портретов при различных видах нелинейности

.
Первое, что надо сделать, это преобразовать данное нелинейное уравнение в дифференциально алгебраическую систему (ДАУ):

.
Это позволяет использовать стандартные методы решения ДАУ: метод Эйлера (неявный) или методы Рунге-Кутты (неявные). Последний более точный и дает возможность оптимизации шага интегрирования.
Эти методы встроены в конкретные алгоритмические языки, например
Python. Подробности я отослал ТС в почте.