clear all; close all;

dx=0.1;
x=-10:dx:10;

y=sech(x);

%exact value of derivative
yd_exact=-sech(x).*tanh(x);

%Plot function and its derivative

figure(1)
set(gca,'FontSize',18)
plot(x,y,'--','LineWidth',2); hold on
plot(x,yd_exact,'LineWidth',2) ;
legend('y(x)','yd exact(x)')

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)
set(gca,'FontSize',18)
plot(x,yd_exact,'o'); hold on
plot(x,yd,'LineWidth',3);
plot(x,yd2)
legend('yd exact(x)','yd order 2','yd order 4')


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

int_trap=trapz(x,y.^2)

int_trap_cum=cumtrapz(x,y.^2)
figure(3)
set(gca,'FontSize',18)
plot(x,int_trap_cum,'LineWidth',2)


sec_sqd_fun=inline('sech(x).^2')

int_quad=quad(sec_sqd_fun,-10,10)