xi = 0 ;
xf = 10 ;
h  = (xf-xi)/1000;
[x,y]=BVPsolve(xi,yesatta(xi),xf,yesatta(xf),10);
plot(x,y,'.-b','LineWidth',1.5);
hold;
xe = [xi:h:xf] ;
plot(xe,yesatta(xe),'-r','LineWidth',0.5);
legend('Numerica','Esatta');