clear all; close all;

dx=0.1;
x=-10:dx:10;
y=sech(x);
yderv=-sech(x).*tanh(x);
plot(x,y,x,yderv)

n=length(x)
% 2nd order accurate
yd(1)=(-3*y(1)+4*y(2)-y(3))/(2*dx);
for j=2:n-1
    yd(j)=(y(j+1)-y(j-1))/(2*dx);
end
yd(n)=(3*y(n)-4*y(n-1)+y(n-2))/(2*dx);

% 4th order accurate
yd2(1)=(-3*y(1)+4*y(2)-y(3))/(2*dx);
yd2(2)=(-3*y(2)+4*y(3)-y(4))/(2*dx);
for j=3:n-2
 yd2(j)=(-y(j+2)+8*y(j+1)-8*y(j-1)+y(j-2))/(12*dx);
end
yd2(n-1)=(3*y(n-1)-4*y(n-2)+y(n-3))/(2*dx);
yd2(n)=(3*y(n)-4*y(n-1)+y(n-2))/(2*dx);


figure(2), plot(x,yderv,'o',x,yd,'c',x,yd2,'m')

% Integration of sech(x)^2
format long e

int_trap=trapz(x,y.^2)

int_quad=quad(inline('sech(x).^2'),-10,10)