Terminal Velocity Differential Equation with from sympy

I am trying to use python's sympy library to calculate from the point of: m(dv/dt) = -mg-kv

I have the following code:

from sympy import *
import matplotlib as plt


mass, g, k, t = symbols('mass g k t', real=True)
v = symbols('v', cls=Function)

# equation = mass*v(t).diff(t)-mass*g+k*(v(t))**2
# print(equation)
equation = Eq(mass*v(t).diff(t), mass*g - k*v(t))
print(equation)

solution = dsolve(equation,v(t),hint='lie_group')
print(solution)

particularSolution = solution.subs([(v(t),0),(t,0)])
print(particularSolution)

C1 = symbols('C1')

const = solve(particularSolution, C1)
print(const)


final_solution = (solution.subs(C1,const[1]).rhs).expand().simplify()
print(final_solution)

I am given the example of a skydiver who's terminal velocity is 200km/hr which I then apply with

#m = 100kg
#g = 9.8 m/s2  =>  127008 km/hr2
#v_max = 200km/hr
#k = mg/v_max   =>   63504 

plot(final_solution.subs([(g,127008),(mass,100),(k,63504)]),(t,0,1))

But the resulting graph looks like

I don't know what is wrong with my equation or sympy implementation to not get an accurate graph

Admittedly, I didn't re-derive your equations before writing this, but at first glance, it appears as though your solutions kinda-sorta converge to what you expect. The only thing that strikes me as strange is that it reaches the limit after one time step. To me, this signals that your initial conditions are not quite correct.

If I put everything back into SI units (meters, kilograms, seconds) and run your code, you seem to get something that smoothly approaches terminal velocity in a realistic amount of time:

mm = 100 #kg
gg = 9.8 #m/s2
vm = 55.5555777778 #m/s
kk = mm*gg/vm 


plot(final_solution.subs([(g,gg),(mass,mm),(k,kk)]),(t,0,100))

which yields

A quick google search of "terminal velocity graphs" show similar results, so I think this is what you are looking for.