Green´s Functions with Reflection
Alberto Cabada, José Ángel Cid,
F. Adrián Fernández Tojo and Beatriz Máquez-Villamarín
Last update: October 2014  on Mathematica 8.0.1.0

Clear["Global`*"];
mess="done";
result="";

CLength[x_]:=Module[{y},y=x;
While[y[[Length[y]]]==0,y=y[[1;;Length[y]-1]]];
Length[y]]

NFill[x_,n_]:=If[TrueQ[n>Length[x]],Join[x,Table[0,{n-Length[x]}]],x];

Start[c1_,c2_,T_,cc1_]:=Module[{asdf,aa,n,bb,bcn,rango2,m,Graphic,opred,Opp,cadenatexto,equation,ecuacion,condcont,condcont2,ecinicial,ecu,eqaux,ViaLibre,c,lc,c1b,c2b},Off[];
ecu=0;
aa=-T;
bb=T;
mess="Processing data...";
n=Max[CLength[c1],CLength[c2]]-1;
If[TrueQ[c1∈Reals&&c2∈Reals&&T∈Reals],Graphic=True,Graphic=False];
ViaLibre=True;
If[Not[n∈Integers&&n>0],MessageDialog["Order must be a positive integer"];
ViaLibre=False;];
If[Not[T>0]&&ViaLibre,MessageDialog["T must be a positive real number"];
ViaLibre=False;];
If[Not[n+1==Length[c1]&&n+1==Length[c2]]&&ViaLibre,MessageDialog["Vector of coefficients: LENGTH INCORRECT"];
ViaLibre=False;];
L[f_][x_]:=Sum[c1b[[k+1]] Derivative[k][f][-x]+c2b[[k+1]] Derivative[k][f][x],{k,0,n}];
R[f_][x_]:=Sum[c1b[[k+1]] Derivative[k][f][-x]-(-1)^k c2b[[k+1]] Derivative[k][f][x],{k,0,n}];
Clear[aux];
Opp=1;
If[Opp==1,aux[var_]:=cc1/.u->var,aux[var_]:=cc1/.{u->var,T->-T}];
EG=False;
If[Not[TrueQ[Norm[c1]*Norm[c2]==0]],m=2 n;
c1b=NFill[c1,m+1];
c2b=NFill[c2,m+1];
bcn=2*CLength[cc1];
c=Table[0,{m+1}];
Do[c[[j+1]]=Sum[(-1)^i*(c1b[[i+1]]*c1b[[j-i+1]]-c2b[[i+1]]*c2b[[j-i+1]]),{i,0,j}],{j,0,m}];
aux2[u_]:=Join[aux[u],Expand[aux[R[u]]]];,m=n;
c1b=NFill[c1,m+1];
c2b=NFill[c2,m+1];
EG=True;
bcn=CLength[cc1];
aux2[var_]:=aux[var];
If[TrueQ[Norm[c1]==0],c=c2;];
If[TrueQ[Norm[c2]==0],c=c1;
Opp=-1;];];
lc=CLength[c]-1;
If[TrueQ[c[[m+1]]==0],MessageDialog["The reduced problem is of order less than 2 n"];
G[t_,s_]="Undetermined";,If[ViaLibre==True,Do[alfa[i,j]=Coefficient[aux2[u][[i]],Derivative[j][u][-T]],{j,0,m-1},{i,1,bcn}];
Do[beta[i,j]=Coefficient[aux2[u][[i]],Derivative[j][u][T]],{j,0,m-1},{i,1,bcn}];
Do[U[i][u_]=Sum[alfa[i,j]*Derivative[j][u][-T]+beta[i,j]*Derivative[j][u][T],{j,0,m-1}],{i,1,bcn}];
condcont2=Sort[Table[Expand[U[i][u]]==0,{i,1,bcn}]];
condcont=Sort[Table[Expand[aux2[u][[i]]]==0,{i,1,bcn}]];
If[TrueQ[Chop[condcont]==Chop[condcont2]],opred[u_][t_]:=Sum[c[[k+1]]Derivative[k][u][t],{k,0,lc}];
D0[k_][u_]:=Derivative[k][u][0];
equation=Join[{opred[y][t]==0},Table[D0[i][y]==0,{i,0,lc-2}],{D0[lc-1][y]==1}];
mess="Solving homogeneous equation...";
ecinicial=DSolve[equation,y,t];
cadenatexto=ToString[ecinicial];
If[Not[StringMatchQ[cadenatexto,"*Root*"]],Result=Style[Column[{Style[" ",Bold],Style["PROBLEM: ",Bold],Style[" ",Bold],Row[{TraditionalForm[L[u][t]]==σ[t],",    t ∈ [",-T,",",T,"]"}],Style[" ",Bold],Style["with boundary conditions",Bold],Style[" ",Bold],Table[aux[u][[i]]==0,{i,1,Length[aux[u]]}],Style[" ",Bold],Style["The Green´s function is giving by: ",Bold],r=ComplexExpand[Re[y/.ecinicial[[1]]]];,mess="Computing fundamental matrix...";
If[TrueQ[c∈Reals&&T∈Reals],Do[soluci[k]=DSolve[Join[{opred[y][t]==0},Table[D0[i][y]==0,{i,0,k-2}],{D0[k-1][y]==1},Table[D0[i][y]==0,{i,k,lc-1}]],y[t],t];
yk[k][t_]=FullSimplify[ComplexExpand[y[t]/.soluci[k][[1]]]];,{k,1,lc}];,Do[yk[k][t_]=Sum[c[[lc+1-j]] Derivative[j-k][r][t],{j,k,lc}];,{k,1,lc}];];
rango2=MatrixRank[Table[U[i][yk[j]],{i,1,bcn},{j,1,lc}]];
If[TrueQ[Not[rango2==lc]],MessageDialog["There is no Green's function for the reduced problem"];
Graphic=False;
G[t_,s_]="There is no unique solution";,eqaux=Table[Sum[beta[i,j]*Derivative[j][r][T-s],{j,0,lc-1}]+Sum[d[j][s] U[i][yk[j]],{j,1,lc}]==0,{i,1,bcn}];
ecuacion=Solve[eqaux,Table[d[j][s],{j,1,lc}]];
If[ecuacion=={},MessageDialog["There is no Green's function"];
Graphic=False;
G[t_,s_]="There is no unique solution";,ecu=1;
mess="Constructing Green's function... (100 s max)";
asdf=ecuacion/.Rule[a_,b_]:>b;
Do[e[j][s_]=d[j][s]/.{d[j][s]->asdf[[1]][[j]]},{j,1,lc}];
h[t_,s_]=Simplify[Sum[e[i][s]*yk[i][t],{i,1,lc}],TimeConstraint->15];
G1[t_,s_]=Simplify[TrigFactor[Chop[r[t-s]]+h[t,s]],TimeConstraint->15];
G2[t_,s_]=Simplify[TrigFactor[h[t,s]],TimeConstraint->15];
Gb[s_][t_]=Piecewise[{{G1[t,s],-T≤s≤T&&-T≤t≤T&&s≤t},{G2[t,s],-T≤s≤T&&-T≤t≤T&&t<s},{0,-T>s||-T>t||T<s||T<t}}];
If[Not[EG],Gb1[t_,s_]=PiecewiseExpand[R[Gb[s]][t],TimeConstraint->15];
G[t_,s_]=Piecewise[{{Simplify[Gb1[t,s],-T≤s≤T∧-T≤t≤T∧s-t≤0∧s+t≤0],-T≤s≤T∧-T≤t≤T∧s-t≤0∧s+t≤0},{Simplify[Gb1[t,s],-T≤s≤T∧-T≤t≤T∧s-t>0∧s+t≤0],-T≤s≤T∧-T≤t≤T∧s-t>0∧s+t≤0},{Simplify[Gb1[t,s],-T≤s≤T∧-T≤t≤T∧s-t≤0∧s+t>0],-T≤s≤T∧-T≤t≤T∧s-t≤0∧s+t>0},{Simplify[Gb1[t,s],-T≤s≤T∧-T≤t≤T∧s-t>0∧s+t>0],-T≤s≤T∧-T≤t≤T∧s-t>0∧s+t>0}}];
G[t_,s_]=Chop[PiecewiseExpand[G[t,s]/c[[m+1]],TimeConstraint->15]];,If[Opp==1,G[t_,s_]=Chop[PiecewiseExpand[Gb[s][t]/c[[m+1]],TimeConstraint->15]];,G[t_,s_]=Chop[PiecewiseExpand[Gb[-s][t]/c[[m+1]],TimeConstraint->15]];];];];];
Row[{Style["G[t,s]= ",Bold],TraditionalForm[G[t,s]]}],Style[" ",Bold],Style[" ",Bold],Style[" ",Bold],If[TrueQ[Graphic],If[ecu==1,Plot3D[G[t,s],{s,aa,bb},{t,aa,bb}],Print["Cannot show the graphic"]]]},{Frame->True,Alignment->Center}]],MessageDialog["Green's Function with a complex expression"];];,MessageDialog["The boundary conditions are not valid"];];];];
mess="done";];

F[]:=(c1={1,0,1};c2={0,0,0};lim=1;cc1={u[1],u[-1]};Nap=False;Framed[Column[{Style["Program to compute the Green's  function of the equation:    ",Bold],Style[Row[{TraditionalForm[Sum[Subscript[a,j] Derivative[j][u][-t],{j,0,n}]+Sum[Subscript[b,j] Derivative[j][u][t],{j,0,n}]==σ[t]],",    t∈[-T,T]"}],Bold],Style[Column[{"                                      "},Center],Bold],Style["with boundary conditions:    ",Bold],Style[Row[{TraditionalForm[Subscript[U,i][u]=Sum[Subsuperscript[α,i,j]Superscript[u,(j)][-T],{j,0,n-1}]+Sum[Subsuperscript[β,i,j]Superscript[u,(j)][T],{j,0,n-1}]==0],",    i=1,...,n"}],Bold],Style[" ",Bold],Column[{Panel[Grid[{{"Coefficients \!\(\*SubscriptBox[\(a\), \(i\)]\)",InputField[Dynamic[c1]]},{"Coefficients \!\(\*SubscriptBox[\(b\), \(i\)]\)",InputField[Dynamic[c2]]},{"T",InputField[Dynamic[lim]]},{"Boundary conditions",InputField[Dynamic[cc1]]},{"Numerical Approximation",Checkbox[Dynamic[Nap],Appearance->Large]}},Alignment->{{Right,Left}}]],Button["Enter",If[Nap,If[Element[c1,Reals],c1=N[c1]];
If[Element[c2,Reals],c2=N[c2]];];Start[c1,c2,lim,cc1],ImageSize->150,Method->"Queued"]},Alignment->Center],Column[{Framed[Style[Row[{"Progress: ",Dynamic[mess]}],Bold]],Dynamic[Result]},Alignment->Left]}]])
F[]