siyyam.m

n=16;
[x,w]=pd(n);t=x2t(n,x);D=deriv(n);X=mult(n);
S=funm(X,@sin);E=expm(X);
A=D^4+S*(D^3+eye(n))+E*(D^2+D);
A(n-3,:)=ones(1,n);
A(n-2,:)=sum(D);A(n-1,:)=sum(D^2);A(n,:)=sum(D^3);
B=t'\exp(-x);
B(n-3)=1;B(n-2)=-1;B(n-1)=1;B(n)=-1;
c=A\B;
rez=(t'*c)-exp(-x);
ind=find(x<10);
plot(x(ind),rez(ind));