h=0.05; x=0:h:pi; y=sin(x)./(x+cos(2*x));
y11=diffctr(y,h,1,1);  y12=diffctr(y,h,1,2);
y21=diffctr(y,h,2,1);  y22=diffctr(y,h,2,2);
y31=diffctr(y,h,3,1);  y32=diffctr(y,h,3,2);
y41=diffctr(y,h,4,1);  y42=diffctr(y,h,4,2);

D=x+cos(2*x); DD=D.^2; DDD=D.^3; D4=D.^4;
E=1-2*sin(2*x); EE=E.^2; EEE=E.^3; E4=E.^4;
yy1=diff(y)/h; f1=cos(x)./D-sin(x).*E./DD;
f2=4*cos(2*x).*sin(x)./DD-y-2*cos(x).*E./DD+2*sin(x).*EE./DDD;
f3=12*cos(x).*cos(2*x)./DD-cos(x)./D-24*cos(2*x).*y.*E./DD+...
   3*y.*E./D+6*cos(x).*EE./DDD-6*y.*EEE./DDD-8*y.*sin(2*x)./D;
f4=96*cos(2*x).^2.*y./DD+y-40*cos(2*x).*y./D+4*cos(x).*E./DD-...
   96*cos(x).*cos(2*x).*E./DDD+4*cos(x).*E./DD+...
   144*cos(2*x).*y.*EE./DDD-12*y.*EE./DDD-24.*cos(x).*EEE./D4+...
   24*y.*E4./D4-32*cos(x).*sin(2*x)./DD+64*y.*E.*sin(2*x)./DD;

L=length(y); i1=3:L; i2=4:L-2; i3=3:L-2; i4=4:L-4;
subplot(221),plot(x,f1,x(i1),y11(1:end-1),'--',x(i2),y12(1:end-1),':')
subplot(222),plot(x,f2,x(i1),y21(1:end-1),'--',x(i2),y22(1:end-1),':')
subplot(223),plot(x,f3,x(i3),y31(1:end-1),'--',x(i4),y32(1:end-1),':')
subplot(224),plot(x,f4,x(i3),y41(1:end-1),'--',x(44),y42(1:end-1),':')