# Plotting a solution of a differential equation with Sagemath

I need to solve a differential equation. The solution will depend on $$t$$ and $$q$$, and I need to define that $$q$$ piecewise depending on $$t$$.

var('k,Tmax,Tmin,w,T0,q'); T=function('T')(t); Te=function('Te')(t);

assume(k>0); assume(Tmax>Tmin); Te(t)=(Tmax+Tmin)/2+(Tmax-Tmin)/2*sin(w*t);


Now this is the differential equation solution:

sol=desolve(diff(T(t),t)-q+k*(T(t)-Te(t)),[T,t],[0,T0]);


The solution with $$q=0$$ for example would be

sol.subs(Tmax=21.6,Tmin=15.2,k=0.024,q=0,T0=15.6,w=pi/12);


but I need that q to be a model for a heater that's on from 8 AM to 22 PM, and off from 22 PM to 8 AM. So I need to define a $$q$$ function that if $$t mod 24$$ is between $$8$$ and $$22$$ its value is $$0.0504$$, and $$0$$ otherwise. Something like this

$$q(t)=\begin{cases}0.0504 \quad\ \ \ \ \ \ \ if \quad t\ mod\ 24 \in[8,22] \\ 0 \ \ \qquad \qquad otherwise\end{cases}$$

I've been trying with the piecewise function but it's not plotting anything, I always get error messages.

You can express $$q(t)$$ as a sum of differences of step functions. Also, it's more efficient to solve the differential equation numerically. I assume you want to plot the solution for some number of days (which can be specified in the code).

days=3
Tmax=21.6; Tmin=15.2; k=0.024; T0=15.6; w=pi/12
var('t');
Te=(Tmax+Tmin)/2+(Tmax-Tmin)/2*sin(w*t)


We define $$q(t)$$:

q = 0.0504*sum(unit_step(t-8-d*24) - unit_step(t-22-d*24) for d in range(days))
plot(q,t,0,days*24)


We define the ODE (also the one with $$q(t)=0$$, to compare):

T=function('T')(t);
ode0 = diff(T,t) == -k*(T-Te)
ode1 = diff(T,t) == q-k*(T-Te)


Finally we solve and plot:

sol0=desolve_rk4(ode0, T, ivar=t, ics=[0,T0], step=0.1, end_points=[0,days*24], output='plot', xmin=0,xmax=days*24)
sol1=desolve_rk4(ode1, T, ivar=t, ics=[0,T0], step=0.1, end_points=[0,days*24], output='plot', xmin=0,xmax=days*24, color='red')
sol0 + sol1