Quit[]

xa=0.0; xb=2.0;
f[x_]:=x*Exp[-4x]
Plot[f[x],{x,xa,xb}]

n=5;
x[1]=xa;
x[n]=xb;
dx=(xb-xa)/(n-1);
Do[x[k]=x[k-1]+dx,{k,2,n}]
Do[a[k]=f[x[k]];h[k]=dx,{k,1,n}]
pts=Table[{x[k],a[k]},{k,1,n}]
Do[alpha[i]=3.*(a[i+1]-a[i])-3.*(a[i]-a[i-1]),{i,2,n-1}];

aa=Table[0.0,{i,1,n},{j,1,n}];
bb=Table[0.0,{j,1,n}];
cc=Table[c[j],{j,1,n}];
MatrixForm[aa]
MatrixForm[bb]
MatrixForm[cc]

(* SPLINE NATURALE *)
Do[aa[[j,j-1]]=dx;aa[[j,j]]=4.0*dx;aa[[j,j+1]]=dx;bb[[j]]=alpha[j]/dx,{j,2,n-1}];
bb[[1]]=0.0;
bb[[n]]=0.0;
aa[[1,1]]=1.0;
aa[[n,n]]=1.0;
MatrixForm[aa]
MatrixForm[bb]

sys=Thread[Dot[aa,cc]==bb];
sol=Solve[sys,cc][[1]];
Do[d[j]=((c[j+1]-c[j])/(3.0*h[j]))/.sol;b[j]=(a[j+1]-a[j])/h[j]-(h[j]*(2.0*c[j]+c[j+1])/3.0)/.sol,{j,1,n-1}]
Do[s[j][y_]=a[j]+b[j]*(y-x[j])+(c[j]/.sol)*(y-x[j])^2+d[j]*(y-x[j])^3,{j,1,n-1}]
Θ[y_]:=If[y>0.0,1.,0.];
fspl[y_]=Sum[s[j][y]*Θ[y-x[j]]*Θ[x[j+1]-y],{j,1,n-1}];
Plot[{fspl[y],f[y]},{y,xa,xb},PlotStyle->{{Thickness[0.005],RGBColor[1,0,0]},{Thickness[0.005],RGBColor[0,0,1]}}]

(* SPLINE COMPLETA *)
Do[aa[[j,j-1]]=dx;aa[[j,j]]=4.0*dx;aa[[j,j+1]]=dx;bb[[j]]=alpha[j]/dx,{j,2,n-1}];
bb[[1]]=3.*(a[2]-a[1])/dx-3.*f'[xa];
bb[[n]]=3.*f'[xb]-3.*(a[n]-a[n-1])/dx;
aa[[1,2]]=dx;
aa[[1,1]]=2.0*dx;
aa[[n,n]]=2.0*dx;
aa[[n,n-1]]=dx;
MatrixForm[aa]
MatrixForm[bb]

sys=Thread[Dot[aa,cc]==bb];
sol=Solve[sys,cc][[1]];
Do[d[j]=((c[j+1]-c[j])/(3.0*dx))/.sol;b[j]=(a[j+1]-a[j])/dx-(dx*(2.0*c[j]+c[j+1])/3.0)/.sol,{j,1,n-1}]
Do[s[j][y_]=a[j]+b[j]*(y-x[j])+(c[j]/.sol)*(y-x[j])^2+d[j]*(y-x[j])^3,{j,1,n-1}]
Θ[y_]:=If[y>0.0,1.,0.];
fspl[y_]=Sum[s[j][y]*Θ[y-x[j]]*Θ[x[j+1]-y],{j,1,n-1}];
Plot[{fspl[y],f[y]},{y,xa,xb}]