# Maple code for Newton Cooling maple lab 2, S2007
#
#Notes on L2.1
 # maple integral table lookup
 unassign('k','omega','t','x');
 integrand:=(35-13*cos(omega*(x-3)))*k*exp(k*x-k*t);
 F:=unapply(integrand,(x,t,k,omega));
 int(F(x,t,k,omega),x=0..t);

 # Test LHS=RHS for u'+ku=kA.
 unassign('t','u0','omega','k'):
 myANS:=your hand-derived formula for u=uh+up:
 LHS:=diff(myANS,t)+k*myANS:
 RHS:=k*(35-13*cos(omega*(t-3))):
 simplify(expand(LHS-RHS));

#Notes on L2.2
 # Extract steady-state solution
 expr:=exp(-k*t)*cos(omega*(t-3))+exp(-2*k*t)*sin(omega*(t-3))
       +2*sin(omega*(t-3));
 SS:=subs(exp(-k*t)=0,exp(-2*k*t)=0,expr);

#Notes on L2.3
 unassign('t','u0','k','omega'):
 AA:=unapply(35-13*cos(omega*(t-3)),(t,omega));
 uss:=35-(13*k/(k^2+omega^2))*(k*cos(omega*(t-3))
      +omega*sin(omega*(t-3))):
 uss0:=subs(t=0,uss);
 U:=unapply((u0-uss0)*exp(-k*t)+uss,(t,u0,k,omega)):

 plot({U(t,74,0.35,Pi/12),AA(t,Pi/12)},t=0..48);

 plot(sin(x),x=0..Pi,color=black);

#Notes on L2.4
with(plots): unassign('t','u0','k','omega'):
 U:=(t,u0,k,omega)->your answer of uh+up from L2.1:
 implicitplot(U(t,74,k,Pi/12)=29,t=0..72,k=0.2..0.45);
 plot3d({U(t,74,k,Pi/12),29},t=0..72,k=0.2..0.45);