Esto sí es amar, matmet

%derivación 5 puntos
syms x
disp('Fórmula de 5 puntos');
f=input('Escriba la función: ');
h=input('Escriba el valor de h: ');
x0=input('Escriba el valor de x0: ');
j=diff(f);
fprintf('Seleccione 1 o 2 según la fórmula que desea usar:nt1- Cinco términos nt2- Cuatro términos (fórmula centrada)')
op=input('nopción: ');
if op==1
        fprintf('n')
        x1=x0+h;
        fprintf('x1= x0+h= %9.15fn',double(x1))
        x2=x0+2*h;
        fprintf('x2= x0+2*h= %9.15fn',double(x2))
        x3=x0+3*h;
        fprintf('x3= x0+3*h= %9.15fn',double(x3))
        x4=x0+4*h;
        fprintf('x4= x0+4*h= %9.15fn',double(x4))
        yaprox=1/(12*h)*(-25*subs(f,x0)+48*subs(f,x1)-36*subs(f,x2)+16*subs(f,x3)-3*subs(f,x4));
        fprintf('Valor aproximado: (1/12h)*(-25fx0+48fx1-36fx2+16fx3-3fx4)) = %9.15fn',double(yaprox))
        yex=subs(j,x0);
        e=abs(yaprox-yex);
end
if op==2
        fprintf('n')
        x1=x0-2*h;
        fprintf('x1= x0-2*h= %9.15fn',double(x1))
        x2=x0-h;
        fprintf('x2= x0-h= %9.15fn',double(x2))
        x3=x0+h;
        fprintf('x3= x0+h= %9.15fn',double(x3))
        x4=x0+2*h;
        fprintf('x4= x0+2*h= %9.15fn',double(x4))
        yaprox=1/(12*h)*(subs(f,x1)-8*subs(f,x2)+8*subs(f,x3)-subs(f,x4));
        fprintf('Valor aproximado = (1/12h)(fx1-8fx2+8fx3-fx4)= %9.15fn',double(yaprox))
        yex=subs(j,x0);
        e=abs(yaprox-yex);
               
end
fprintf('El valor exacto es: %9.15fn',double(yex));
fprintf('El error es: %9.15fn',double(e));

%Derivación 3 puntos________________________________________________________________________________________________________________
%derivación // Hecho por Diego Quintanilla
syms x
disp('Fórmula de 3 puntos');
f=input('Escriba la función: ');
h=input('Escriba el valor de h: ');
x0=input('Escriba el valor de x0: ');
j=diff(f);
fprintf('Seleccione 1 o 2 según la fórmula que desea usar:nt1-Tres puntos con tres términos nt2-Tres puntos con dos términos (fórmula centrada)')
op=input('nopción: ');
if op==1
        fprintf('n')
        x1=x0+h;
        fprintf('x1= x0+h = %9.15fn',double(x1));
        x2=x0+2*h;
        fprintf('x2= x0+2*h = %9.15fn',double(x2));
        yaprox=1/(2*h)*(-3*subs(f,x0)+4*subs(f,x1)-subs(f,x2));
        fprintf('Valor aproximado= (1/2h)(-3fx0+4fx1-fx2) = %9.15fn',double(yaprox))
        yex=subs(j,x0);
        e=abs(yaprox-yex);
end
if op==2
        fprintf('n')
        x1=x0+h;
        fprintf('x1=x0+h = %9.15fn',double(x1));
        x2=x0-h;
        fprintf('x1=x0-h = %9.15fn',double(x2));
        yaprox=1/(2*h)*(subs(f,x1)-subs(f,x2));
        fprintf('El valor aproximado es: (1/2h)(fx1-fx2) = %9.15fn',double(yaprox));
        yex=subs(j,x0);
        e=abs(yaprox-yex);
               
end
fprintf('El valor exacto es: %9.15fn',double(yex));
fprintf('El error es: %9.15fn',double(e));

%________________________________________________________________________________________________________________________________

%funcion entre comillas, o punto y coma al final

   disp('1-->Extrapolacion con funcion');
   disp('2-->Extrapolacion con tabla');
   select = input('seleccion : ');
   arreglo = zeros();
  switch(select)
       case 1
          syms x
           f=input('digite la funcion f(x) : ' );
           N=input('digite la cantidad de iteraciones N(j) : ');
           h=input('digite h : ');
           x0=input('digite el punto inicial x0 : ');
           fprintf('nn');
           for i=1: N
               div=h;
               d=1;
               if(i~=1)
                 div = (h/2^(i-1));
                 d=2^(i-1);
               end
               superior = x0+div;
               inferior = x0-div;
               arreglo(1,i) = (1/(2*div))*(subs(f,superior)-subs(f,inferior));
               fprintf('N(%1.0f)(h/%1.0f)=(1/(2h/%1.0f))*(F(x0+h/%1.0f)-F(x0-h/%1.0f))=%3.15fnn',1,d,d,d,d,arreglo(1,i));
           end
           for i=1:N
               m=1;
               fprintf('n---------------------------------------------------------------------------------nn');
               for j=1:(N-i)
                   k=2^(j);
                   d=(4^(i)-1);
                   arreglo((i+1) , j)=(arreglo(i,(j+1)))+((arreglo(i,(j+1))-arreglo(i,j))/d);
                   fprintf('nn N(%1.0f)(h/%1.0f)=N(%1.0f)(h/%1.0f)+((N(%1.0f)(h/%1.0f)-N(%1.0f)(h/%1.0f))/%1.0f =%3.15f',(i+1),m,i,k,i,k,(j+1),m,d,arreglo((i+1),j));
                   m=m*2;
                   Vaprox = arreglo((i+1),j);
               end
           end
           syms x;
           dif=diff(f,x);
           Vreal = double(subs(dif , x0));
           fprintf('nValor aproximado es =%9.15f ',Vaprox);
           fprintf('nValor exacto es =%9.15f ' , Vreal );
           err=abs(Vreal-Vaprox);
           fprintf('nError =%3.3e nn',err);
       case 2
           
       otherwise
           return;
  end
%n(2,2)=n(2,1)+[n(2,1)-n(1,1)]/k   
%n(3,2)=n(3,1)+[n(3,1)-n(2,1)]/k 
%n(4,2)=n(4,1)+[n(4,1)-n(3,1)]/k

%____________________________________________________________________________________________________________________________________

%Integracion de Romberg
format long
syms x  y  z;
g=input('Escriba el integrando g(x)=');
n=input('Escriba el tamaño R(n,n) n: ');
a=input('limite a=');
b=input('limite b=');
RM=zeros(n,n);
fa=subs(g,a);
fb=subs(g,b);
for j=1:n
    h(j)=(b-a)/(2^(j-1));
    fprintf('Valores para h cuando j=1:n')
    h(j)
end
RM(1,1)=(h(1)/2)*(fa+fb);
for j=2:n
    Sum=0;
    for i=1:(2^(j-2))
        c=(a+((2*i-1)*h(j)));
        fc=subs(g,c);
        Sum=Sum+fc;
    end
    RM(j,1)=(1/2)*(RM(j-1,1)+h(j-1)*Sum);
end
for j=2:n
    for i=j:n
        RM(i,j)=RM(i,j-1) + (RM(i,j-1)-RM(i-1,j-1))/((4^(j-1))-1);
    end
end
for j=1:n
    for i=j:n
        if j==1&&i==1
            fprintf('R(1,1)=h1/2[fa+fb]=%9.15fn',RM(1,1))
        elseif j==1&&i>1
            fprintf('R(%1.0f,1)=(1/2)[R(%1.0f,1)+h%1.0f SumF(a+(2i-1)h%1.0f]=%9.15fn',i,i-1,i-1,i,RM(i,1))
        else
        fprintf('R(%1.0f,%1.0f)=R(%1.0f,%1.0f)+R(%1.0f,%1.0f)-R(%1.0f,%1.0f)/(4^(%1.0f-1)-1)=%9.15fn',i,j,i,j-1,i,j-1,i-1,j-1,j,RM(i,j))
        end
    end
end
fxd=input('Escriba la misma función pero en términos de y, f(y): ');
iex=int(fxd,y,a,b);
e=abs(RM(n,n)-iex);
fprintf('nValor aproximado: %9.15fn',RM(n,n))
fprintf('nValor real: %9.15fn',double(iex))
fprintf('nEl error es: %9.15fn',double(e))

%_________________________________________________________________________________________________________________________
function D=nca(a,b,n) %Newton Cotes abierta
disp('Newton-Cotes abiertas')
syms x
F=input('Escriba la función f(x): ');
f=inline(F);
h=(b-a)/(n+2);
fprintf('Seleccione la fórmula a aplicar: nt 1-Regla del punto medio (n=0) nt 2-Para n=1 nt 3-Para n=2 nt 4-Para n=3n');
op=input('opción: ');

if op==1
    x0=a+h;
    fx0=f(x0);
    inte=2*h*fx0;
    iex=int(F,x,a,b);
    e=abs(iex-inte);
     
end
if op==2
    x0=a+h;
    x1=x0+h;
    fx0=f(x0);
    fx1=f(x1);
    inte=3*h/2*(fx0+fx1);
    iex=int(F,x,a,b);
    e=abs(iex-inte);
        
end
if op==3
    x0=a+h;
    x1=x0+h;
    x2=x0+2*h;
    fx0=f(x0);
    fx1=f(x1);
    fx2=f(x2);
    inte=4*h/3*(2*fx0-fx1+2*fx2);
    iex=int(F,x,a,b);
    e=abs(iex-inte);
    
end
if op==4
    x0=a+h;
    x1=x0+h;
    x2=x0+2*h;
    x3=x0+3*h;
    fx0=f(x0);
    fx1=f(x1);
    fx2=f(x2);
    fx3=f(x3);
    inte=5*h/24*(11*fx0+fx1+fx2+11*fx3);
    iex=int(F,x,a,b);
    e=abs(iex-inte);
       
end
    fprintf('nEl valor aproximado es: %9.15f',double(inte));
    fprintf('nEl valor exacto es: %9.15f',double(iex));
    fprintf('nEl error es: %9.15fn',double(e));
end
%______________________________________________________________________________________________________________________
function D=ncc(a,b,n)% Newton Cotes cerradas
disp('Newton-Cotes cerradas')
syms x
F=input('Escriba la función f(x): ');
f=inline(F);
h=(b-a)/n;
fprintf('Seleccione la fórmula a aplicar: nt 1-Regla del trapecio nt 2-Regla de Simpson nt 3-Regla de 3/8 de Simpson nt 4-Regla de Boolen');
op=input('opción: ');

if op==1
    fa=f(a);
    fb=f(b);
    inte=h/2*(fa+fb);
    iex=int(F,x,a,b);
    e=abs(iex-inte);
     
end
if op==2
    fa=f(a);
    fb=f(b);
    x1=a+h;
    fx1=f(x1);
    inte=h/3*(fa+4*fx1+fb);
    iex=int(F,x,a,b);
    e=abs(iex-inte);
        
end
if op==3
    fa=f(a);
    fb=f(b);
    x1=a+h;
    x2=a+2*h;
    fx1=f(x1);
    fx2=f(x2);
    inte=3*h/8*(fa+3*fx1+3*fx2+fb);
    iex=int(F,x,a,b);
    e=abs(iex-inte);
    
end
if op==4
    fa=f(a);
    fb=f(b);
    x1=a+h;
    x2=a+2*h;
    x3=a+3*h;
    fx1=f(x1);
    fx2=f(x2);
    fx3=f(x3);
    inte=2*h/45*(7*fa+32*fx1+12*fx2+32*fx3+7*fb);
    iex=int(F,x,a,b);
    e=abs(iex-inte);
    
    
end
    fprintf('nEl valor aproximado es: %9.15f',double(inte));
    fprintf('nEl valor exacto es: %9.15f',double(iex));
    fprintf('nEl error es: %9.15fn',double(e));
end
%____________________________________________________________________________________________________________________________________

%Compuesta de Simpson 
fprintf('>>clcn');
fprintf('>>syms x;n');
n=input('>>n=');
a=input('>>a=');
b=input('>>b=');
disp('>>h=(b-a)/n;');
h=(b-a)/n;
disp('>>x0=a');
a
for i=1:n
   fprintf('>>x%1.0f=(a+%1.0f*h)n',i,i);
   z=a+(i)*h;
   fprintf('n x%1.0f=nnt%1.15fnn',i,z);
end
disp('>>x12=b');
b
syms x;
g=input('>>g=');
iex=int(g,x,a,b);
g
sum1=0;
sum2=0;
cad1='';
cad2='';
for i=1:n-1
    y=a+i*h;
    if(mod(i,2)==0)
        sum1 = sum1 + subs(g,y);
        nod=strcat('subs(g,x' ,num2str(i) ,')'); 
        cad1 = strcat(cad1 , '+' , nod);
    else
        nod2=strcat('subs(g,x' ,num2str(i) ,')'); 
        sum2= sum2 + subs(g,y);
        cad2 = strcat(cad2 , '+' , nod2);
    end
end
fprintf('>>sum1=%s n',cad1);
sum1
fprintf('n');
fprintf('>>sum2=%s n',cad2);
sum2
fprintf('n');
fa=subs(g,a);
disp('>>fa=subs(g,a);');
fb=subs(g,b);
disp('>>fb=subs(g,b);');
fprintf('>>faprox=(h/3)*(fa + 2*sum1 + 4*sum2 + fb)n');
faprox=(h/3)* (fa + 2*sum1 + 4*sum2 + fb);
e=abs(iex-faprox);
fprintf('El valor aproximado es: %9.15fn',double(faprox));
fprintf('El valor exacto es: %9.15fn',double(iex));
fprintf('El error es: %9.15fn',double(e));

%___________________________________________________________________________________________________________

%Trapecio compuesto:
function D=trapeciocomp(a,b,n)
syms y
f=input('Escriba la función f(y): ');
fa=subs(f,a);
fb=subs(f,b);
h=(b-a)/n;
iex=int(f,y,a,b);
%función suma = fs
fs=0;
    for i=1:n-1;
        x(i)=a+h*i;
        fs=fs+subs(f,x(i));
        
    end
    
    inte=0.5*h*(fa+2*fs+fb);
    e=abs(inte-iex);
    fprintf('El valor aproximado de la integral es: %9.15fn',double(inte));
    fprintf('El valor exacto de la integral es: %9.15fn',double(iex));
    fprintf('El error es: %9.15fn',double(e));
end

%_________________________________________________________________________________________________________________________
%Cuadratura Gaussiana
format long;
P=input('>> Escriba Polinomio como vector [1 0 ...]=');
N=input('>> Escriba cuánto vale N=: ');
disp('>>r=roots(P)');
r=roots(P); display(r); c=zeros(N);
syms x;
for i=1:N
    fprintf('>>x%1.0f=%9.15f n', i, r(i));
end
for i=1:N
   num=1;
   denom=1;
   cad='';
   fprintf('>>c%1.0f=',i);
   for j=1:N
       if(i~=j)
           num= num * (x-r(j));
           denom= denom * (r(i)-r(j));
           a=['x',num2str(i)];
           b=['x',num2str(j)];
           if(j~= N)
                cad=strcat(cad ,strcat('(',a, '-' ,b , ')' , '*'));
                fprintf('(x-x%1.0f)*', j);
           else
                fprintf('(x-x%1.0f)', j);
                cad=strcat(cad ,strcat('(',a, '-' ,b , ');'));
           end
       end
   end
   fprintf('/(%s)', cad);
   fprintf('n>>int(c%1.0f, -1 , 1 , x )n',i);
   c(i)=int((num/denom),x,-1,1);
   c(i)
end
syms x t
%Respuestas del polinomio
f=input('Escriba el integrando: ');
pol=input('Escriba el polinomio f(x): ');
a=input('a= ');
b=input('b= ');
T=double(solve(pol,x));
n=length(T);
format long
%Integracion de polinomios de Lagrange
for i=1:n
    M=1;
    for j=1:n
        if(i==j)
        else
            M=M*((x-T(j))/(T(i)-T(j)));
        end
    end
    C(i)=double(int(M,x,-1,1));
end
%Sustitucion de la funcion integrando a la nueva funcion
V=(1/2)*(((b-a)*t)+a+b);
dV=(b-a)/2;
ft=subs(f,V)*dV;

S=0;
%Uso de las C y de la nueva funcion evaluada en las raices del polinomio
for i=1:n
   S=S+(C(i)*subs(ft,T(i))); 
end
double(S)
%Valor exacto
Exacto=(int(f,x,a,b));
fprintf('El valor aproximado es %9.15f n',double(S));
fprintf('El valor exacto es %9.15f n',double(Exacto));
fprintf('El error %e n',double(abs(S-Exacto)));

%___________________________________________________________________________________________________________________________________

%Adaptivo de cuadratura
function D=adapcua(a,b)
syms x
f=input('Escriba la función f(x): ');
h=(b-a)/2;
fa=subs(f,a);
fb=subs(f,b);
fx1=subs(f,a+h/2);
fx2=subs(f,a+h);
fx3=subs(f,a+3*h/2);
iaprox=(h/6)*(fa+4*fx1+fx2)+(h/6)*(fx2+4*fx3+fb);
fprintf('El valor aproximado de la integral es: %9.15fn',double(iaprox))
end
%__________________________________________________________________________________________________________________________________

%Segunda derivada
disp('Seguda derivada')
f=input('Escriba la función: ');
syms x
h=input('Escriba el valor de h: ');
x0=input('Escriba el valor de x0: ');
sd=1/h^2*(subs(f,x0-h)-2*subs(f,x0)+subs(f,x0+h));
y=diff(f);
yy=diff(y);
y2=subs(yy,x0);
fprintf('El valor aproximado de la segunda derivada es: %9.15fn',double(sd))
fprintf('El valor exacto de la segunda derivada es: %9.15fn',double(y2));

%____________________________________________________________________________________________________________________________

%Compuesta de Simpson para integrales impropias singularidad límite
%izquierdo, cuando es derecho, hacer cambio x=z;
%Ahorita está diseñado para n=16;
function D=compimpr(a,b,n,p,sl)
h=(b-a)/n;
H=double(h)
x0=a;
%Aquí agregar o quitar x según sea el tamaño de n
x1=x0+h;x2=x1+h;x3=x2+h;x4=x3+h;x5=x4+h;x6=x5+h;x7=x6+h;x8=x7+h;
x9=x8+h;x10=x9+h;x11=x10+h;x12=x11+h;x13=x12+h;x14=x13+h;x15=x14+h;x16=x15+h;

syms x;
g=input('Función g(x): ');
d1=diff(g);d2=diff(d1);d3=diff(d2);d4=diff(d3);
g1=subs(g,sl); g2=subs(d1,sl); g3=subs(d2,sl); g4=subs(d3,sl); g5=subs(d4,sl);
px=g1+g2*(x-(sl))+g3*(x-(sl))^2/2+g4*(x-(sl))^3/6+g5*(x-(sl))^4/24;
f=(g-px)/(x-(sl))^p;

%Aquí agregar los valores de x
sum1=subs(f,x2)+subs(f,x4)+subs(f,x6)+subs(f,x8)+subs(f,x10)+subs(f,x12)+subs(f,x14);%suma los pares - el último
sum2=subs(f,x1)+subs(f,x3)+subs(f,x5)+subs(f,x7)+subs(f,x9)+subs(f,x11)+subs(f,x13)+subs(f,x15);%Suma todo impar

sumatoria1=double(sum1)
sumatoria2=double(sum2)
w=px/(x-(sl))^p;
%Sí el límite con sigularidad es el a (izquierdo), no evaluar subs(f,a) dentro de i1.
%Si es el b (derecho), no evaluar subs(f,b);
i1=double(h/3*(2*sum1+4*sum2+subs(f,b)))
i2=double(int(w,x,a,b))
iaprox=double(i1+i2);
fprintf('El valor aproximado es: %9.15fn',double(iaprox))

end
%____________________________________________________________________________________________________________________________

%Derivadas
%función    derivada
%asin(u)=u'/sqrt(1-u^2)
%acos(u)=-u'/sqrt(1-u^2)
%atan(u)=u'/(1+u^2)
%acot(u)-u'/(1+u^2)
%asec(u)=u'/(abs(u)*sqrt(u^2  -1)
%acsc(u)=-u'/(abs(u)*sqrt(u^2  -1)
%a^u=a^u(ln(u))u'
%sinx = cosx
%cosx = -sinx
%tanx = sec^2
%cotx = -csc^2
%secx = secx tanx
%cscx = -cscx cotx

%integración
%sqrt(a^2-u^2)            sqrt(a^2+u^2)                    sqrt(u^2-a^2)
%u=a sin            u=a tan                     u= a sec
du= a cos                       du=a sec^2                                      du= a tan sec
sqrt = a cos                    sqrt= a sec                                     sqrt= a tan

%1/sqrt(a^2-u^2) = arcosin(u/a)
%1/(a^2 +u^2)=(1/a)(arcotan(u/a))
%1/(u*sqrt(u^2 - a^2) = 1/a(arcosec(u/a))

%longitud de arco__________________________________
%fórmula=sqrt(1+^[f'(x)]^2)

%area de superficie de revolución_____________________
%formula= 2pi*f(x)*sqrt(1+^[f'(x)]^2)

%trabajo________________________________
%formula= (densidad)(volumen pequeño diferencial)(distancia) 
%nota: bombeo es llenar, vaciar es distancia disco a bomba (k-y) llena es Y

%fuerza________________________________________
%formula= densidad*h(y)L(y)

%Arandela________________________________
%formula= pi(R^2  -   r^2)
%eje a la izquierda o abajo                 eje derecha o arriba
%f(x) - eje                                   eje -f(x)

%Capa cilindrica_____________________________
%formula= 2pi*p(x)h(x) nota: p(x)=x; h(x)= funcion de tope rec,tipico
%eje a la izquierda o abajo                 eje derecha o arriba
%f(x) - eje                                   eje -f(x)

%método de discos_____________________________--
%formula= pi[f(x)]^2