This question has already been answered. So disregard this question. I don't know how to delete the question.

function traction_ode45_a clf, clear all global grade m den Cd Af g fo fi eff PT m=2500; den=1.218; Cd=0.35; Af-1.5; g=9.81; f0-0.0136; f1=5.184e-7; eff-0.8; Vi=0.1; grade=4; PT-60000; [t60, V60]=ode45 (@velocity, [0 1501, Vi); grade=4; PT-80000; [t80, 180] Fode45 (@Velocity, [0 150), Vi); grade=4; PT-100000; [t100, V100)-ode 45 (@Velocity, [0 150), Vi); grade=4; PT-120000; [t120, V120)-ode 45 (@Velocity, [0 150), Vi); rt-3600/1000; plot (t60, V60*rt, t80,V80*rt, t100, V100*rt, t120, V120*rt), grid xlabel('Time, sec'), ylabel('Velocity, km/hr') title(['m', num2str (m), 'kg, Cd', num2str(Cd),', Af-', num2str (AC), ' m2, PT-', num2str (PT/1000), 'kW']) legend('60 kW', '80 kW', '100 kW', '120 kW', 'location', 'northwest') function vdot -Velocity(t,v) global grade m den Cd Af g fo fl eff PT theta-atan (grade/100); W-m*g; r1-0.5*den+CdAf+f1*W*cos (theta): r2-W* (f0*cos (theta)+sin(theta)); r3=PT+eff; vdot-(1v.^3+r2*v-r3)./(m*v)