syms g m b F
syms v(t)
eq=m*diff(v(t),t)+b*v(t)==F

vel=simplify(dsolve(eq,v(0)==0))

%

u = symunit;
unitConversionFactor(u.N, u.kg*u.m/u.s^2);
% b*v(t) on Njuutonites, leiame b ühiku
syms b_units
b_units = simplify(solve(b_units*(u.m/u.s)==u.N))

%
% Paneme parameetrid sisse ja arvutame lõppkiiruse

vel_SI = subs(vel,[m,b,F],[1000*u.kg, 50*b_units, 1400*u.N])

%
% Lõppkiirus piirväärtusega ja m/s ning ft/s avaldisena

vel_final=limit(vel_SI,t,Inf), 
vel_final=rewrite(vel_final,u.m)
vel_final_Imperial=rewrite(vel_final,u.ft)

%
% Plottimiseks määrame ajaühiku sekundites ja eraldame valemi ning ühikud

syms T, vel_s=subs(vel_SI,t,T*u.s);
[vel_formula, units_SI]=separateUnits(vel_s);
units_SI=rewrite(units_SI,u.m);
fplot(vel_formula,[0 100])
xlabel('Time in s')
ylabel(['Velocity in ' symunit2str(units_SI)])