######################################################################
##Math421: Save this file as Math421. To use it, stay in the         #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read Math421 : <Enter>                             #
##(in Windows (DOS) you need to type:                                #
##read `C:\\direct1\\direct2\\Math421`:                              #
#assuming that Math421 resides in C:\direct1\direct2
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg@math.rutgers.edu.                                          # 
#######################################################################

#Created: Nov. 25, 2002
#This version: Dec. 5, 2002 (More Heat Eq. stuff and ODE)
#Previous version: Dec. 2, 2002 (Thanks to Philip Engel!)
#Math421: A Maple package to accopmany Math421
#Please report bugs to zeilberg@math.rutgers.edu

print(`To save your work in dumb terminals use the command`):
print(`writeto; Use the Maple help to find out about it`):
print(``):
print(`Created: Nov. 25, 2002`):
print(`Version of Dec. 5, 2002`):
print(`Previous version: Dec. 2, 2002 `):
print(`(Thanks to Philip Engel, for correcting)`):
print(` numerous bugs and erros!`):
print(`This version of Dec. 5 also has several kinds of  Heat Eqs`):
print(`And Stability and Nature of Critical Points (sections 10.3, 10.4)`):
print(`The exercises refer to the textbook, Advanced Engineering Mathematics`):
print(`by Peter V. O'Neil, fifth edition `):
lprint(``):
print(`Written by Doron Zeilberger, zeilberg@math.rutgers.edu`):
lprint(``):
print(`Please report bugs to zeilberg@math.rutgers.edu`):
lprint(``):
 print(`The most current version of this  package `):
 print(` is available from`):
 print(`http://www.math.rutgers.edu/~zeilberg`):
 print(`For a list of the procedures type Help(), for help with`):
 print(`a specific procedure, type Help(procedure_name)`):

print(`Math421, a Maple package to accompany Math421 Dr. Z.'s class`):
print(`This is for MANDATORY self-study computer-labs`):
print(`To replace Thu. Nov. 19`):
print(`MANDATORY HOMEWORK: For each of the procedures listed`):
print(`by typing  Help();  find TWO  problems`):
print(`from ALL the homework assigned in this class`):
print(`that can be answered by that procedure.`):

Help:=proc()
if args=NULL then
 print(`Contains the following procedures: Critical,`):
 print(` dAlembert, FourierSeries `):
 print(` FourierSineSeries, FourierCosineSeries , `):
 print(` HeatEqDifferentTemp, HeatEqInsulatedEnds, `):
 print(`InHomogWaveEq,  WaveEq, WaveEqZeroDisp, WaveEqZeroInitVel `):
fi:

if nops([args])=1 and op(1,[args])=Critical then
print(`Critical(a,b,c,d): Given a 2 by 2 matrix A=[[a,b],[c,d]]`):
print(`finds the nature and stability and asympt. stability of`):
print(`the system x'(t)=a*x(t)+b*y(t), y'(t)=c*x(t)+d*y(t)`):
print(`For example, to do ex. 10.3.6 (and 10.4.6) do`):
print(`Critical(2,-7,5,-10):`):
fi:

if nops([args])=1 and op(1,[args])=dAlembert then
print(`dAlembert(f,g,x,t,c): d'Alembert's solution of the Wave Eq`):
print(`u_tt=c^2 u_xx for -infinity<x<infinity, t>0`):
print(`with u(x,0)=f, u_t(x,0)=g(x) -infinity<x<infinity, `):
print(`For example, to do Ex. 16.4.2 type`):
print(` dAlembert(x^2-2*x,cos(x),x,t,4);`):
fi:

if nops([args])=1 and op(1,[args])=ExtractSines then
print(`ExtractSines(f,x,L): extracts the pure building blocks`):
fi:

if nops([args])=1 and op(1,[args])=ExtractSinesWithExp then
print(`ExtractSinesWithExp(f,x,t,k,L): `):
print(`Given an expression with sin(n*x*Pi/L) in it`):
print(`extracts those with exp(-k*(n*Pi/L)^2*t) attached, followed by`):
print(` the rest without the sines `):
fi:

if nops([args])=1 and op(1,[args])=FourierSeries then
print(`FourierSeries(f,x,L): The Fourier series of a function`):
print(`of x given by the expression f, over the interval [-L,L] `):
print(`For example try FourierSeries(x^2,x,Pi);`):
fi:

if nops([args])=1 and op(1,[args])=FourierCosineSeries then
print(`FourierCosineSeries(f,x,L): The Fourier Cosine series of a function`):
print(`of x given by the expression f, over the interval [0,L] `):
print(`For example try FourierCosineSeries(x^2,x,Pi);`):
fi:

if nops([args])=1 and op(1,[args])=FourierSineSeries then
print(`FourierSineSeries(f,x,L): The Fourier Sine series of a function`):
print(`of x given by the expression f, over the interval [0,L] `):
print(`For example try FourierSineSeries(x^2,x,Pi);`):
fi:

if nops([args])=1 and op(1,[args])=HeatEqDifferentTemp then
print(`HeatEqDifferentTemp(f,k,x,t,L,T1,T2): The solution of the Heat `):
print(` equation`):
print(`u_t=k u_xx (0<x<L,t>0) with u(0,t)=T1, u(L,t)=T2 (t>=0), `):
print(` u(x,0)=f(x) for 0<=x<=L. For example to do `):
print( ` Ex. 17.2.17, type : HeatEqDifferentTemp(x^2,16,x,t,1,2,5); `):
fi:

if nops([args])=1 and op(1,[args])=HeatEqInsulatedEnds then
print(`HeatEqInsulatedEnds(f,k,x,t,L): The solution of the Heat equation`):
print(`u_t=k u_xx (0<x<L,t>0) with u_x(0,t)=u_x(L,t)=0 (t>=0), `):
print(`u(x,0)=f(x) for 0<=x<=L. For example to do`):
print(`Ex. 17.2.4, type : HeatEqInsulatedEnds(sin(x),1,x,t,Pi);`):
fi:

if nops([args])=1 and op(1,[args])=HeatEqZeroTemp then
print(`HeatEqZeroTemp(f,k,x,t,L): The solution of the Heat equation`):
print(`u_t=k u_xx (0<x<L,t>0) with u(0,t)=u(L,t)=0 (t>=0), `):
print(`u(x,0)=f(x) for 0<=x<=L. For example to do`):
print(`Ex. 17.2.1, type : HeatEqZeroTemp(x*(L-x),k,x,t,L);`):
fi:

if nops([args])=1 and op(1,[args])=InHomogWaveEq then
print(`InHomogWave(F,f,g,x,t,c): Solution of`):
print(`u_tt(x,t)=c^2 u_xx (x,t)+ F(x,t) for -infinity<x<infinity, t>0`):
print(`with u(x,0)=f, u_t(x,0)=g(x) -infinity<x<infinity, `):
print(`For example, to do Ex. 16.4.20 type`):
print(` InHomogWaveEq(t*sin(x),1-x^2,x*sin(x),x,t,c);`):
fi:

if nops([args])=1 and op(1,[args])=WaveEq then
print(`WaveEq(f,g,c,x,t,L): The solution of the wave equation`):
print(`y_tt=c^2 y_xx in [0,L] with y(0,t)=y(L,t)=0, y(x,0)=f and`):
print(`y_t(x,0)=g (0<x<L). For example to do`):
print(`Ex. 16.2.9, type :`):
print(`WaveEq(x^2*(x-4),1,sqrt(8),x,t,4);`):
fi:

if nops([args])=1 and op(1,[args])=WaveEqZeroDisp then
print(`WaveEqZeroDisp(g,c,x,t,L): The solution of the wave equation`):
print(`y_tt=c^2 y_xx in [0,L] with y(0,t)=y(L,t)=0, y(x,0)=0 and`):
print(`y_t(x,0)=g (0<x<L). For example to do`):
print(`Ex. 16.2.1, type :`):
print(`WaveEqZeroDisp(2*x*(Heaviside(x)-Heaviside(x-1)),c,x,t,2);`):
fi:

if nops([args])=1 and op(1,[args])=WaveEqZeroInitVel then
print(`WaveEqZeroInitVel(f,c,x,t,L): The solution of the wave equation`):
print(`y_tt=c^2 y_xx in [0,L] with y(0,t)=y(L,t)=0, y(x,0)=f(x) and`):
print(`y_t(x,0)=0 (0<x<L). For example to do`):
print(`Ex. 16.2.2, type :WaveEqZeroInitVel(2*sin(Pi*x),3,x,t,4);`):
fi:

end:

Simp:=proc(A,n)
local A1,A2,sol,i,gu:
A2:=denom(A):
A1:=numer(A):
sol:={solve(A2,n)}:
gu:=[]:
for i from 1 to nops(sol) do
 if type(sol[i],integer) and sol[i]>0 then
  gu:=[op(gu),[sol[i],limit(A,n=sol[i])]]:
 fi:
od:
gu:

end:

#FourierSeries(f,x,L): The Fourier series of a function
#of x given by the expression f, over the interval [-L,L] 
#For example try FourierSeries(x^2,x,Pi);
FourierSeries:=proc(f,x,L)
local a0,an,bn,n,F1,g,f1,F2:

F1:=ExtractSines(f,x,L):
F2:=ExtractCosines(f,x,L):
g:=F1[1]+F2[1]:
f1:=expand(f-g):

assume(n,integer):
a0:=(1/L)*simplify(int(f1,x=-L..L)):
an:=simplify((1/L)*simplify(int(f1*cos(n*\Pi*x/L),x=-L..L))):
bn:=simplify((1/L)*simplify(int(f1*sin(n*\Pi*x/L),x=-L..L))):
if an=0 and bn=0 then
RETURN(a0/2+g):
elif an<>0 and bn=0 then
RETURN(a0/2+g+Sum(an*cos(n*Pi*x/L),n=1..infinity)):
elif an=0 and bn<>0 then
RETURN(g+1/2*a0+Sum(bn*sin(n*Pi*x/L),n=1..infinity)):
else
RETURN(g+1/2*a0+
Sum(an*cos(n*Pi*x/L),n=1..infinity)+Sum(bn*sin(n*Pi*x/L),n=1..infinity)):
fi:



end:

#FourierSineSeries(f,x,L): The Fourier Sine series of a function
#of x given by the expression f, over the interval [0,L] 
#For example try FourierSineSeries(x^2,x,Pi);
FourierSineSeries:=proc(f,x,L)
local bn,n,F,g,f1:
F:=ExtractSines(f,x,L):
g:=F[1]:
f1:=F[2]:
assume(n,integer):
bn:=2/L*int(f1*sin(n*\Pi*x/L),x=0..L):
bn:=simplify(bn):
if bn=0 then
RETURN(g):
else
RETURN(g+Sum(bn*sin(n*Pi*x/L),n=1..infinity)):
fi:
end:

#FourierCosineSeries(f,x,L): The Fourier Cosine series of a function
#of x given by the expression f, over the interval [0,L] 
#For example try FourierCosineSeries(x^2,x,Pi);
FourierCosineSeries:=proc(f,x,L)
local a0, an,n,F, g,f1:
F:=ExtractCosines(f,x,L):
g:=F[1]:
f1:=F[2]:

a0:=2/L*int(f1,x=0..L):
assume(n,integer):
an:=2/L*int(f1*cos(n*\Pi*x/L),x=0..L):
an:=simplify(an):
if an<>0 then
RETURN(a0/2+g+Sum(an*cos(n*Pi*x/L),n=1..infinity)):
else
RETURN(a0/2+g):
fi:
end:


#ExtractSines(f,x,L): Given an expression with sin(n*x*Pi/L) in it
#extracts those, followed by the rest
ExtractSines:=proc(f,x,L)
local i,f1,f11,f111,g,c,F:
f1:=f:
F:=f1:
g:=0:

if type(f1,`+`) then
for i from 1 to nops(f1) do
f11:=op(i,f1):

if type(f11,`*`) then
 c:=op(1,f11):
 f111:=op(2,f11):
   if type(c,numeric) and type(f111,function) then
      if op(0,f111)=sin and type(op(1,f111)/(Pi*x/L),integer) then
         g:=g+f11:
         F:=F-f11:
      fi:
   fi:
fi:

if type(f11,function) then
  c:=1:
 f111:=f11:
   if type(c,numeric) and type(f111,function) then
      if op(0,f111)=sin and type(op(1,f111)/(Pi*x/L),integer) then
         g:=g+f11:
         F:=F-f11:
      fi:
   fi:
fi:

od:
fi:

if not type(f1,`+`) then
f11:=f1:

if type(f11,`*`) then
 c:=op(1,f11):
 f111:=op(2,f11):
   if type(c,numeric) and type(f111,function) then
      if op(0,f111)=sin and type(op(1,f111)/(Pi*x/L),integer) then
         g:=g+f11:
         F:=F-f11:
      fi:
   fi:
fi:

if type(f11,function) then
  c:=1:
 f111:=f11:
   if type(c,numeric) and type(f111,function) then
      if op(0,f111)=sin and type(op(1,f111)/(Pi*x/L),integer) then
         g:=g+f11:
         F:=F-f11:
      fi:
   fi:
fi:


fi:

g,F:
end:


#ExtractCosines(f,x,L): Given an expression with cos(n*x*Pi/L) in it
#extracts those, followed by the rest
ExtractCosines:=proc(f,x,L)
local i,f1,f11,f111,g,c,F:
f1:=f:
F:=f1:
g:=0:

if type(f1,`+`) then
for i from 1 to nops(f1) do
f11:=op(i,f1):

if type(f11,`*`) then
 c:=op(1,f11):
 f111:=op(2,f11):
   if type(c,numeric) and type(f111,function) then
      if op(0,f111)=cos and type(op(1,f111)/(Pi*x/L),integer) then
         g:=g+f11:
         F:=F-f11:
      fi:
   fi:
fi:

if type(f11,function) then
  c:=1:
 f111:=f11:
   if type(c,numeric) and type(f111,function) then
      if op(0,f111)=cos and type(op(1,f111)/(Pi*x/L),integer) then
         g:=g+f11:
         F:=F-f11:
      fi:
   fi:
fi:

od:
fi:

if not type(f1,`+`) then
f11:=f1:

if type(f11,`*`) then
 c:=op(1,f11):
 f111:=op(2,f11):
   if type(c,numeric) and type(f111,function) then
      if op(0,f111)=cos and type(op(1,f111)/(Pi*x/L),integer) then
         g:=g+f11:
         F:=F-f11:
      fi:
   fi:
fi:

if type(f11,function) then
  c:=1:
 f111:=f11:
   if type(c,numeric) and type(f111,function) then
      if op(0,f111)=cos and type(op(1,f111)/(Pi*x/L),integer) then
         g:=g+f11:
         F:=F-f11:
      fi:
   fi:
fi:


fi:

g,F:
end:

         
            
#ExtractSinesWithCosines(f,x,t,c,L): 
#Given an expression with sin(n*x*Pi/L) in it
#extracts those with cosines attached, followed by the rest
#(without the cosines)
ExtractSinesWithCosines:=proc(f,x,t,c,L)
local i,f1,f11,f111,g,c1,F,mumu:
f1:=f:
F:=f1:
g:=0:

if type(f1,`+`) then
for i from 1 to nops(f1) do
f11:=op(i,f1):

if type(f11,`*`) then
 c1:=op(1,f11):
 f111:=op(2,f11):
   if type(c1,numeric) and type(f111,function) then
    mumu:=op(1,f111)/(Pi*x/L) :
      if op(0,f111)=sin and type(mumu,integer) then
         g:=g+f11*cos(mumu*Pi*c*t/L):
         F:=F-f11:
      fi:
   fi:
fi:

if type(f11,function) then
  c1:=1:
 f111:=f11:
    mumu:=op(1,f111)/(Pi*x/L) :
   if type(c1,numeric) and type(f111,function) then
      if op(0,f111)=sin and type(mumu,integer) then
         g:=g+f11*cos(mumu*Pi*c*t/L):
         F:=F-f11:
      fi:
   fi:
fi:

od:
fi:

if not type(f1,`+`) then
f11:=f1:

if type(f11,`*`) then
 c1:=op(1,f11):
 f111:=op(2,f11):
   if type(c1,numeric) and type(f111,function) then
     mumu:=op(1,f111)/(Pi*x/L):
      if op(0,f111)=sin and type(mumu,integer) then
         g:=g+f11*cos(mumu*Pi*c*t/L):
         F:=F-f11:
      fi:
   fi:
fi:

if type(f11,function) then
  c1:=1:
 f111:=f11:
   if type(c1,numeric) and type(f111,function) then
            mumu:=op(1,f111)/(Pi*x/L):
      if op(0,f111)=sin and type(mumu,integer) then
         g:=g+f11*cos(mumu*Pi*c*t/L):
         F:=F-f11:
      fi:
   fi:
fi:


fi:

g,F:
end:

#WaveEqZeroInitVel(f,c,x,t,L): The solution of the wave equation
#y_tt=c^2 y_xx in [0,L] with y(0,t)=y(L,t)=0, y(x,0)=f(x) and
#y_t(x,0)=0 (0<x<L). For example to do
#Ex. 16.2.2, type :WaveEqZeroInitVel(2*sin(Pi*x),3,x,t,4);
WaveEqZeroInitVel:=proc(f,c,x,t,L)
local bn,n,F,g,f1:
F:=ExtractSinesWithCosines(f,x,t,c,L):
g:=F[1]:
f1:=F[2]:
assume(n,integer):
bn:=2/L*int(f1*sin(n*\Pi*x/L),x=0..L):
bn:=simplify(bn):
if bn=0 then
RETURN(g):
else
RETURN(g+Sum(bn*sin(n*Pi*x/L)*cos(n*Pi*c*t/L),n=1..infinity)):
fi:
end:



#ExtractSinesWithSines(f,x,t,c,L): 
#Given an expression with sin(n*x*Pi/L) in it
#extracts those with sines attached, followed by the rest
#(without the cosines)
ExtractSinesWithSines:=proc(f,x,t,c,L)
local i,f1,f11,f111,g,c1,F,mumu:
f1:=f:
F:=f1:
g:=0:

if type(f1,`+`) then
for i from 1 to nops(f1) do
f11:=op(i,f1):

if type(f11,`*`) then
 c1:=op(1,f11):
 f111:=op(2,f11):
   if type(c1,numeric) and type(f111,function) then
    mumu:=op(1,f111)/(Pi*x/L) :
      if op(0,f111)=sin and type(mumu,integer) then
         g:=g+f11*sin(mumu*Pi*c*t/L):
         F:=F-f11:
      fi:
   fi:
fi:

if type(f11,function) then
  c1:=1:
 f111:=f11:
    mumu:=op(1,f111)/(Pi*x/L) :
   if type(c1,numeric) and type(f111,function) then
      if op(0,f111)=sin and type(mumu,integer) then
         g:=g+f11*sin(mumu*Pi*c*t/L):
         F:=F-f11:
      fi:
   fi:
fi:

od:
fi:

if not type(f1,`+`) then
f11:=f1:

if type(f11,`*`) then
 c1:=op(1,f11):
 f111:=op(2,f11):
   if type(c1,numeric) and type(f111,function) then
     mumu:=op(1,f111)/(Pi*x/L):
      if op(0,f111)=sin and type(mumu,integer) then
         g:=g+f11*sin(mumu*Pi*c*t/L):
         F:=F-f11:
      fi:
   fi:
fi:

if type(f11,function) then
  c1:=1:
 f111:=f11:
   if type(c1,numeric) and type(f111,function) then
            mumu:=op(1,f111)/(Pi*x/L):
      if op(0,f111)=sin and type(mumu,integer) then
         g:=g+f11*sin(mumu*Pi*c*t/L):
         F:=F-f11:
      fi:
   fi:
fi:


fi:

g,F:
end:

#WaveEqZeroDisp(g,c,x,t,L): The solution of the wave equation
#y_tt=c^2 y_xx in [0,L] with y(0,t)=y(L,t)=0, y(x,0)=0 and
#y_t(x,0)=g (0<x<L). For example to do
#Ex. 16.2.1, type :
#WaveEqZeroDisp(2*x*(Heaviside(x)-Heaviside(x-1)),c,x,t,2);
WaveEqZeroDisp:=proc(g,c,x,t,L)
local bn,n,F,g1,f1:
F:=ExtractSinesWithSines(g,x,t,c,L):
g1:=F[1]:
f1:=F[2]:
assume(n,integer):
bn:=2/(n*Pi*c)*int(f1*sin(n*\Pi*x/L),x=0..L):
bn:=simplify(bn):
if bn=0 then
RETURN(g1):
else
RETURN(g1+Sum(bn*sin(n*Pi*x/L)*sin(n*Pi*c*t/L),n=1..infinity)):
fi:
end:


#WaveEq(f,g,c,x,t,L): The solution of the wave equation
#y_tt=c^2 y_xx in [0,L] with y(0,t)=y(L,t)=0, y(x,0)=f and
#y_t(x,0)=g (0<x<L). For example to do
#Ex. 16.2.9, type :
#WaveEq:=proc(x^2*(x-4),1,sqrt(8),x,t,4);
WaveEq:=proc(f,g,c,x,t,L):
WaveEqZeroInitVel(f,c,x,t,L)+WaveEqZeroDisp(g,c,x,t,L):
end:


#dAlembert(f,g,x,t,c): d'Alembert's solution of the Wave Eq
#u_tt=c^2 u_xx for -infinity<x<infinity, t>0
#with u(x,0)=f, u_t(x,0)=g(x) -infinity<x<infinity, 
#For example, to do Ex. 16.4.2 type
# dAlembert(x^2-2*x,cos(x),x,t,4);
dAlembert:=proc(f,g,x,t,c)
local xi:
(subs(x=x-c*t,f)+subs(x=x+c*t,f))/2+1/2/c*int(subs(x=xi,g),xi=x-c*t..x+c*t):
end:


#InHomogWaveEq(F,f,g,x,t,c): Solution of
#u_tt(x,t)=c^2 u_xx (x,t)+ F(x,t) for -infinity<x<infinity, t>0
#with u(x,0)=f, u_t(x,0)=g(x) -infinity<x<infinity, 
#For example, to do Ex. 16.4.20 type
# InHomogWaveEq(t*sin(x),1-x^2,x*sin(x),x,t,c);
InHomogWaveEq:=proc(F,f,g,x,t,c)
local xi,eta:
(subs(x=x-c*t,f)+subs(x=x+c*t,f))/2+1/2/c*int(subs(x=xi,g),xi=x-c*t..x+c*t)+
1/(2*c)*int(int(subs({x=xi,t=eta},F),xi=x-c*t+c*eta..x+c*t-c*eta),eta=0..t):
end:


#ExtractSinesWithExp(f,x,t,k,L): 
#Given an expression with sin(n*x*Pi/L) in it
#extracts those with exp(-k*(n*Pi/L)^2*t) attached, followed by the rest
#(without the sines)
ExtractSinesWithExp:=proc(f,x,t,k,L)
local i,f1,f11,f111,g,c1,F,mumu:
f1:=f:
F:=f1:
g:=0:

if type(f1,`+`) then
for i from 1 to nops(f1) do
f11:=op(i,f1):

if type(f11,`*`) then
 c1:=op(1,f11):
 f111:=op(2,f11):
   if type(c1,numeric) and type(f111,function) then
    mumu:=op(1,f111)/(Pi*x/L) :
      if op(0,f111)=sin and type(mumu,integer) then
         g:=g+f11*exp(-(mumu*Pi/L)^2*k*t):
         F:=F-f11:
      fi:
   fi:
fi:

if type(f11,function) then
  c1:=1:
 f111:=f11:
    mumu:=op(1,f111)/(Pi*x/L) :
   if type(c1,numeric) and type(f111,function) then
      if op(0,f111)=sin and type(mumu,integer) then
         g:=g+f11*exp(-(mumu*Pi/L)^2*k*t):
         F:=F-f11:
      fi:
   fi:
fi:

od:
fi:

if not type(f1,`+`) then
f11:=f1:

if type(f11,`*`) then
 c1:=op(1,f11):
 f111:=op(2,f11):
   if type(c1,numeric) and type(f111,function) then
     mumu:=op(1,f111)/(Pi*x/L):
      if op(0,f111)=sin and type(mumu,integer) then
         g:=g+f11*exp(-(mumu*Pi/L)^2*k*t):
         F:=F-f11:
      fi:
   fi:
fi:

if type(f11,function) then
  c1:=1:
 f111:=f11:
   if type(c1,numeric) and type(f111,function) then
            mumu:=op(1,f111)/(Pi*x/L):
      if op(0,f111)=sin and type(mumu,integer) then
         g:=g+f11*exp(-(mumu*Pi/L)^2*k*t):
         F:=F-f11:
      fi:
   fi:
fi:


fi:

g,F:
end:

#HeatEqZeroTemp(f,k,x,t,L): The solution of the Heat equation
#u_t=k u_xx (0<x<L,t>0) with u(0,t)=u(L,t)=0 (t>=0), 
#u(x,0)=f(x) for 0<=x<=L. For example to do
#Ex. 17.2.1, type : HeatEqZeroTemp(x*(L-x),k,x,t,L);
HeatEqZeroTemp:=proc(f,k,x,t,L)
local cn,n,F,g1,f1:
F:=ExtractSinesWithExp(f,x,t,k,L):
g1:=F[1]:
f1:=F[2]:
assume(n,integer):
cn:=(2/L)*int(f1*sin(n*\Pi*x/L),x=0..L):
cn:=simplify(cn):
if cn=0 then
RETURN(g1):
else
RETURN(g1+Sum(cn*sin(n*Pi*x/L)*exp(-(n*Pi/L)^2*k*t),n=1..infinity)):
fi:
end:

#ExtractCosinesWithExp(f,x,t,k,L): 
#Given an expression with cos(n*x*Pi/L) in it
#extracts those with exp(-k*(n*Pi/L)^2*t) attached, followed by the rest
#(without the sines)
ExtractCosinesWithExp:=proc(f,x,t,k,L)
local i,f1,f11,f111,g,c1,F,mumu:
f1:=f:
F:=f1:
g:=0:

if type(f1,`+`) then
for i from 1 to nops(f1) do
f11:=op(i,f1):

if type(f11,`*`) then
 c1:=op(1,f11):
 f111:=op(2,f11):
   if type(c1,numeric) and type(f111,function) then
    mumu:=op(1,f111)/(Pi*x/L) :
      if op(0,f111)=cos and type(mumu,integer) then
         g:=g+f11*exp(-(mumu*Pi/L)^2*k*t):
         F:=F-f11:
      fi:
   fi:
fi:

if type(f11,function) then
  c1:=1:
 f111:=f11:
    mumu:=op(1,f111)/(Pi*x/L) :
   if type(c1,numeric) and type(f111,function) then
      if op(0,f111)=cos and type(mumu,integer) then
         g:=g+f11*exp(-(mumu*Pi/L)^2*k*t):
         F:=F-f11:
      fi:
   fi:
fi:

od:
fi:

if not type(f1,`+`) then
f11:=f1:

if type(f11,`*`) then
 c1:=op(1,f11):
 f111:=op(2,f11):
   if type(c1,numeric) and type(f111,function) then
     mumu:=op(1,f111)/(Pi*x/L):
      if op(0,f111)=cos and type(mumu,integer) then
         g:=g+f11*exp(-(mumu*Pi/L)^2*k*t):
         F:=F-f11:
      fi:
   fi:
fi:

if type(f11,function) then
  c1:=1:
 f111:=f11:
   if type(c1,numeric) and type(f111,function) then
            mumu:=op(1,f111)/(Pi*x/L):
      if op(0,f111)=cos and type(mumu,integer) then
         g:=g+f11*exp(-(mumu*Pi/L)^2*k*t):
         F:=F-f11:
      fi:
   fi:
fi:


fi:

g,F:
end:

#HeatEqInsulatedEnds(f,k,x,t,L): The solution of the Heat equation
#u_t=k u_xx (0<x<L,t>0) with u_x(0,t)=u_x(L,t)=0 (t>=0), 
#u(x,0)=f(x) for 0<=x<=L. For example to do
#Ex. 17.2.4, type : HeatEqInsulatedEnds(sin(x),1,x,t,Pi);
HeatEqInsulatedEnds:=proc(f,k,x,t,L)
local cn,n,F,g1,f1,c0:
F:=ExtractCosinesWithExp(f,x,t,k,L):
g1:=F[1]:
f1:=F[2]:
c0:=(2/L)*int(f1,x=0..L):
assume(n,integer):
cn:=(2/L)*int(f1*cos(n*\Pi*x/L),x=0..L):
cn:=simplify(cn):
if cn=0 then
RETURN(g1):
else
RETURN(g1+1/2*c0+Sum(cn*cos(n*Pi*x/L)*exp(-(n*Pi/L)^2*k*t),n=1..infinity)):
fi:
end:


#HeatEqDifferentTemp(f,k,x,t,L,T1,T2): The solution of the Heat equation
#u_t=k u_xx (0<x<L,t>0) with u(0,t)=T1, u(L,t)=T2 (t>=0), 
#u(x,0)=f(x) for 0<=x<=L. For example to do
#Ex. 17.2.17, type : HeatEqDifferentTemp(x^2,16,x,t,1,2,5);
HeatEqDifferentTemp:=proc(f,k,x,t,L,T1,T2)
local psi,f1:
psi:=(T2-T1)/L*x+T1:
f1:=f-psi:
HeatEqZeroTemp(f1,k,x,t,L)+psi:
end:


#Critical(a,b,c,d): Given a 2 by 2 matrix A=[[a,b],[c,d]]
#finds the nature and stability and asympt. stability of
#the system x'(t)=a*x(t)+b*y(t), y'(t)=c*x(t)+d*y(t)
#For example, to do ex. 10.3.6 (and 10.4.6) do
#Critical(2,-7,5,-10):
Critical:=proc(a,b,c,d)
local lam,eq,gu:
eq:=expand((a-lam)*(d-lam)-b*c):
gu:=evalf(evalc({solve(eq,lam)})):

if nops(gu)=1 then
lam:=evalf(gu[1]):
if a=lam and d=lam and b=0 and c=0 then

 if lam<0 then
  RETURN([proper_node,sink, stable,asymptotically_stable]):
  elif lam>0 then
   RETURN([proper_node,source, unstable]):
  fi:

else

 if lam<0 then
  RETURN([improper_node,sink, stable,asymptotically_stable]):
  elif lam>0 then
   RETURN([improper_node,source, unstable]):
  fi:

 
fi:
fi:

if nops(gu)=2 then
 if coeff(gu[1],I,1)=0 and evalf(gu[1]*gu[2])>0 then
    if gu[1]<0 then
     RETURN([nodal_sink,stable, asymptotically_stable]):
    else
     RETURN([nodal_source,unstable]):
    fi:
 elif coeff(gu[1],I,1)=0 and evalf(gu[1]*gu[2])<0 then
       RETURN([saddle_point,unstable]):
 elif coeff(gu[1],I,1)<>0 and evalf(gu[1]-I*coeff(gu[1],I,1))>0 then
       RETURN([spiral_source,unstable]):
 elif coeff(gu[1],I,1)<>0 and evalf(gu[1]-I*coeff(gu[1],I,1))<0 then
        RETURN([spiral_sink,stable,asymptotically_stable]):
 elif coeff(gu[1],I,1)<>0 and evalf(gu[1]-I*coeff(gu[1],I,1))=0 then
        RETURN([center,stable,not_asymptotically_stable]):
 else
   ERROR(`Something is wrong`):
fi:
fi:
end: