(* Ficheiro do Mathematica para a solucao do problema 6 *)

(* Input de variaveis *)
Clear[lbd,eta]

lbd=20
eta=2
c=-3.
d=3.

(* 
 Valor para comecar a encontrar a raiz com FindRoot. E preciso confirmar 
 com os  graficos para ver se faz sentido. O Mathematica encontra as 
 solucoes para y > yini
*)

yini=Pi/2 + 0.001   

(* Fim do input *)


(* Definicao das funcoes para a solucoes par, impar, eta >>1, e Delta y *)

f = Function[{y,lbd}, 2 y Sqrt[lbd-y^2]/(2 y^2 -lbd)]

feven = Function[{y, lbd, eta}, y Sqrt[lbd - y^2]*(1 + 
            Tanh[eta Sqrt[lbd - y^2]])/(y^2 - (lbd - y^2)*
              Tanh[eta*Sqrt[lbd - y^2]])]

fodd = Function[{y, lbd, eta}, y Sqrt[lbd - y^2]*(1 + 
            1/Tanh[eta Sqrt[lbd - y^2]])/(y^2 - (lbd - y^2)/
              Tanh[eta*Sqrt[lbd - y^2]])]

Deltay = Function[{y, lbd, eta}, 4 E^(-2*eta*Sqrt[lbd - y^2])*
      y*(lbd - y^2)/lbd/(2 + Sqrt[lbd - y^2])]


(* Calculo das solucoes: y0, yeven, yodd *)

r0 = FindRoot[f[y, lbd] == Tan[y], {y, yini}]; y0 = y /. r0

reven = FindRoot[feven[y, lbd, eta] == Tan[y], {y, yini}]; yeven = y /. reven

rodd = FindRoot[fodd[y, lbd, eta] == Tan[y], {y, yini}]; yodd = y /. rodd ;

Deltayexact= yodd-yeven;

(* Grafico para observar as solucoes *)

plotodd=Plot[{Tan[y],fodd[y,lbd,eta]},{y,0,Sqrt[lbd]},PlotRange->{c,d}];
ploteven=Plot[{Tan[y],feven[y,lbd,eta]},{y,0,Sqrt[lbd]},PlotRange->{c,d}]

(* Calculo das funcoes de onda par e impar *)

alfa=Function[{y,lbd},Sqrt[lbd-y^2]]

C4=1
C1even=C4
C1odd=-C4
A2even=C4 E^(-(eta+1)*alfa[yeven,lbd])*(Cos[(eta+1)*yeven] +
		alfa[yeven,lbd]/yeven*Sin[(eta+1)*yeven])
B2even=C4 E^(-(eta+1)*alfa[yeven,lbd])*(Sin[(eta+1)*yeven] -
		alfa[yeven,lbd]/yeven*Cos[(eta+1)*yeven])
A1even=A2even
B1even=-B2even
C2even=(A2even*Cos[yeven*eta]+B2even*Sin[yeven*eta])/Cosh[eta*alfa[yeven,lbd]]
A2odd=C4 E^(-(eta+1)*alfa[yodd,lbd])*(Cos[(eta+1)*yodd] +
		alfa[yodd,lbd]/yodd*Sin[(eta+1)*yodd])
B2odd=C4 E^(-(eta+1)*alfa[yodd,lbd])*(Sin[(eta+1)*yodd] -
		alfa[yodd,lbd]/yodd*Cos[(eta+1)*yodd])
A1odd=-A2odd
B1odd=B2odd
C3odd=(A2odd*Cos[yodd*eta]+B2odd*Sin[yodd*eta])/Sinh[eta*alfa[yodd,lbd]]

(* Funcoes nao normalizadas *)

psievenaux=Function[{x,lbd,eta},
If[x<-(eta+1),C1even E^(alfa[yeven,lbd] x),
   If[x<-eta,A1even*Cos[yeven x]+B1even*Sin[yeven x],
      If[x<eta,C2even*Cosh[alfa[yeven,lbd] x],
	 If[x<eta+1,A2even*Cos[yeven x]+B2even*Sin[yeven x],C4 E^(-alfa[yeven,lbd] x)]]]]];

psioddaux=Function[{x,lbd,eta},
If[x<-(eta+1),C1odd E^(alfa[yodd,lbd] x),
   If[x<-eta,A1odd*Cos[yodd x]+B1odd*Sin[yodd x],
      If[x<eta,C3odd*Sinh[alfa[yodd,lbd] x],
	 If[x<eta+1,A2odd*Cos[yodd x]+B2odd*Sin[yodd x],C4 E^(-alfa[yodd,lbd] x)]]]]];

(* Normalizacao *)

Neven=Integrate[psievenaux[x,lbd,eta]^2,{x,-10*eta,10*eta}]
Nodd=Integrate[psioddaux[x,lbd,eta]^2,{x,-10*eta,10*eta}]

(* Funcoes Normalizadas *)

psieven=Function[x,psievenaux[x,lbd,eta]/Sqrt[Neven]]
psiodd=Function[x,psioddaux[x,lbd,eta]/Sqrt[Nodd]];

