%esercizio 2 5.9.22
clear all
close all
clc

x=linspace(0,1);

f=@(x) sin(x) +x.^2 +7*x;
g=@(x) 5*exp(-x);

h=@(x) f(x)-g(x); %fplot(h,[0,1]) tra 0 ed 1 non vale mai zero (stabilità per newton raphson);
dh=@(x) cos(x) + 2*x +5*exp(-x) +7 ;
Nmax=100;
toll=sqrt(eps);

[x0] = NEWTONRAPHSON(h,dh,x,Nmax,toll);
% fplot(h,'r')
% hold on
% plot(x0,h(x0),'*b')

xN=linspace(0,3,5);
yN=h(xN);
[a] = diffdivise(xN,yN);
xx=linspace(0,3);
 
PN1=@(xx) PN(a,xN,xx);
p=polyfit(xN,yN,4);
P=@(xx) polyval(p,xx);

c=0;
d=3;
for i=0:4
    x1(i+1)= (c+d)/2 + (d-c)/2*cos((2*i+1)/(10)*pi);
end

y1=h(x1);

P2=@(xx) Lagrange(x1,y1,xx);

figure
fplot(h,[0,3])
hold on
% fplot(P,[0,3])
fplot(PN1,[0,3])
fplot(P2,[0,3])
plot(xN,yN,'*r')
legend('h','P.Newton','P.CHebyshev')

function [x0] = NEWTONRAPHSON(h,dh,x,Nmax,toll)
%cerca zero
a=x(1);
b=x(end);

x0 = a - sign(h(a))*(b-a)/(sign(h(b))-sign(h(a)));
k=1;
e(k)=abs(h(x0));

while k<Nmax && e(k)>toll
    x0=x0 - h(x0)/dh(x0);
    k=k+1;
    e(k)=abs(h(x0));
end
end

function [a] = diffdivise(xi,yi)
%algoritmo differenze divise
n=length(xi)-1;
a=yi;
for i=1:n
    for j=n:-1:i
        a(j+1)=(a(j+1)-a(j))./(xi(j+1)-xi(j));
    end
end
end


function [y] = PN(a,xi,xx)
%valutazione polinomio di newton
n=length(a);
y=0;

for i=n-1:-1:0
    y= a(i+1) +(xx-xi(i+1)).*y;
end

end

function [h] = Lagrange(xi,yi,xx)
%Lagrange + chebyshev
x=xi;
y=yi;
n=length(xi)-1;
P=xx-xx;

for j=0:n-1
    num=1;
    den=1;
    for i=0:n
        num=num.*(xx-x(i+1));
        den=den.*(x(j+1)-x(i+1));
    end
    P=P+y(j+1).*num./den;
end
h=P;
end