(* PROBLEMI AL CONTORNO *)

(* ESERCIZIO 11.3.1 *)

Clear[a,b,p,q,r,w,x,aa]
fp[x_]:=0;
fq[x_]:=4;
fr[x_]:=-4x;
a=0.;
b=1.;
b=4.;
α=0.;
β=2.;
xrange={a,b+0.2};
eqdiff=y''[x]==fp[x]*y'[x]+fq[x]*y[x]+fr[x];
soldiff=DSolve[{eqdiff,y[a]==α,y[b]==β},y[x],x];
ysol[x_]=y[x]/.soldiff

n=2;
h=(b-a)/n;
x[1]=a;
x[n+1]=b;
ww=Table[w[i],{i,1,n+1}];
Do[x[i+1]=x[i]+h,{i,1,n}];
Do[p[i]=fp[x[i]];q[i]=fq[x[i]];r[i]=fr[x[i]],{i,1,n+1}];
aa=Table[0,{i,1,n+1},{j,1,n+1}];
aa[[1,1]]=1.;
aa[[n+1,n+1]]=1.;
MatrixForm[aa];
Do[aa[[i,i-1]]=2.+h*p[i];aa[[i,i+1]]=2.-h*p[i];aa[[i,i]]=-2.*(2.+h^2*q[i]),{i,2,n}];
MatrixForm[aa]
bb=Table[2.*h^2*r[i],{i,1,n+1}];
bb[[1]]=α;
bb[[n+1]]=β;
MatrixForm[bb]

sol=Solve[Dot[aa,ww]==bb,ww][[1]];
fsol=Table[{x[i],w[i]/.sol},{i,1,n+1}];
pl1=ListPlot[fsol,PlotStyle->{PointSize[0.02],RGBColor[1,0,0]}];
pl2=Plot[ysol[x],{x,a,b}];
Show[pl1,pl2,PlotRange->{xrange,{-0.2,3.2}}]

(* ESERCIZIO 11.3.2 *)

Clear[a,b,p,q,r,w,x,aa]
fp[x_]:=1;
fq[x_]:=2;
fr[x_]:=Cos[x];
a=0.;
b=0.5*Pi;
α=-0.3;
β=-0.1;
eqdiff=y''[x]==fp[x]*y'[x]+fq[x]*y[x]+fr[x];
soldiff=DSolve[{eqdiff,y[a]==α,y[b]==β},y[x],x];
ysol[x_]=y[x]/.soldiff

n=8;
h=(b-a)/n;
x[1]=a;
x[n+1]=b;
ww=Table[w[i],{i,1,n+1}];
Do[x[i+1]=x[i]+h,{i,1,n}];
Do[p[i]=fp[x[i]];q[i]=fq[x[i]];r[i]=fr[x[i]],{i,1,n+1}];
aa=Table[0,{i,1,n+1},{j,1,n+1}];
aa[[1,1]]=1.;
aa[[n+1,n+1]]=1.;
MatrixForm[aa];
Do[aa[[i,i-1]]=2.+h*p[i];aa[[i,i+1]]=2.-h*p[i];aa[[i,i]]=-2.*(2.+h^2*q[i]),{i,2,n}];
MatrixForm[aa]
bb=Table[2.*h^2*r[i],{i,1,n+1}];
bb[[1]]=α;
bb[[n+1]]=β;
MatrixForm[bb]

sol=Solve[Dot[aa,ww]==bb,ww][[1]];
fsol=Table[{x[i],w[i]/.sol},{i,1,n+1}];
pl1=ListPlot[fsol,PlotStyle->{PointSize[0.02],RGBColor[1,0,0]}];
pl2=Plot[ysol[x],{x,a,b}];
Show[pl1,pl2,PlotRange->{-0.4,-0.1}]

(* ESERCIZIO LIBERO *)

Clear[a,b,p,q,r,w,x,aa]
fp[x_]:=0;
fq[x_]:=-9;
fr[x_]:=0;
a=0.;
b=1.;
α=1.;
β=0.5;
xrange={a,b};
yrange={-20,20};
eqdiff=y''[x]==fp[x]*y'[x]+fq[x]*y[x]+fr[x];
soldiff=DSolve[{eqdiff,y[a]==α,y[b]==β},y[x],x];
ysol[x_]=y[x]/.soldiff

n=4;
h=(b-a)/n;
x[1]=a;
x[n+1]=b;
ww=Table[w[i],{i,1,n+1}];
Do[x[i+1]=x[i]+h,{i,1,n}];
Do[p[i]=fp[x[i]];q[i]=fq[x[i]];r[i]=fr[x[i]],{i,1,n+1}];
aa=Table[0,{i,1,n+1},{j,1,n+1}];
aa[[1,1]]=1.;
aa[[n+1,n+1]]=1.;
MatrixForm[aa];
Do[aa[[i,i-1]]=2.+h*p[i];aa[[i,i+1]]=2.-h*p[i];aa[[i,i]]=-2.*(2.+h^2*q[i]),{i,2,n}];
MatrixForm[aa]
bb=Table[2.*h^2*r[i],{i,1,n+1}];
bb[[1]]=α;
bb[[n+1]]=β;
MatrixForm[bb]

sol=Solve[Dot[aa,ww]==bb,ww][[1]];
fsol=Table[{x[i],w[i]/.sol},{i,1,n+1}];
pl1=ListPlot[fsol,PlotStyle->{PointSize[0.02],RGBColor[1,0,0]}];
pl2=Plot[ysol[x],{x,a,b}];
Show[pl1,pl2,PlotRange->{xrange,yrange}]