Hello, I've made a script to solve a non-linear ordinary differential equation for the temperature (lumped-mass) of a brake rotor due to some constant deceleration of a race car (Qbrake in script). The problem I'm having is properly simulating the velocity of the car. What I'd like for it to do is start from 60mph (Vmax), decelerate (@ Gdec=-1.6 g's) to 30mph (Vmin), and accelerate (@ Gacc=0.5 g's) back to 60 mph, and repeat this process for the time span (a few minutes) of ODE45 which should drive the temperature to a roughly steady value as this velocity cycles between Vmin and Vmax. I've attempted this w/ IF conditionals and a persistent "gees" (the acceleration value changing based on velocity conditions), but running the script results in erratic velocity values over the time span. That is, it does not run smoothly between the velocity limits. Would would be a better way to simulate this behavior within the ODE45 algorithm? My script: Code (Text): function thermalrotor % Initial conditions Vo=60*1.67; % mph --> ft/s To=75; % degF, ambient % ODE45 tspan=0:.1:120; [t,y]=ode45(@thermalt,tspan,[Vo To]) Vc=y(:,1)/1.67 Tr=y(:,2); Tm=2795; % melting temperature of ductile iron plot(t,Tr,t,Tm,'g') title('Temperature as a function of time, ductile iron rotor'); ylabel('Temperature, F'); xlabel('Time, s'); legend('Rotor Temp.','Melting Temp.'); clear %% Subfunction function dx=thermalt(t,x) % Constants: persistent gees W=575; % lbf g=32.18; % ft/s Vmax=60*1.67; % ft/s Vmin=30*1.67; % ft/s Gacc=.5; Gdec=-1.6; As1=48.027/144; % ft^2 As2=38.799/144; % ft^2 Ah=4*0.2281/144; % ft^2 Lh=0.31846/12; % ft M=0.9335; % lbm Cp=0.128; % Btu/(lbm*degF) k=46/3600; % Btu/(s*ft*degF) h=25/3600; % Btu/(s*ft^2*degF) fv=1.0; fe=0.9; sigma=1.714*10^-9/3600; % Btu/(ft^2*s*R^4) Tamb=75; % degF rad=fe*fv*sigma*As2; if x(1)>=Vmax gees=Gdec; elseif x(1)<=Vmin gees=Gacc; end % Energy absorbed by rotor from deceleration (Q=Fv): Qbrake=-W/3*x(1)*gees/778; % Energy radiated to air Qrad=rad*((Tamb+460)^4-(x(2)+460)^4); % x(1) = Vc = car velocity (linear) % x(2) = Tr = lumped mass rotor temperature dx=zeros(2,1); dx(1)=gees*g; dx(2)=(Qbrake+(h*As2+k*Ah/Lh)*(Tamb-x(2))+Qrad)/(M*Cp); Thanks!