METODO DE HEUN

%MÉTODO DE HEUN

% - Introduzca la ecuación diferencial        : 'Dy=y-(x^2)+1'
% - Introduzca la condición y(a)=b            : 'y(0)=0.5'
% - Introduzca la función de trabajo          : y-(x^2)+1
% - Introduzca la condición inicial           : 0.5
% - Introduzca el valor de a                  : 0
% - Introduzca el valor de b                  : 1
% - Introduzca el tamaño de paso h            : 0.1

fprintf('\n');
clear all
clc
fprintf('                                     --------------\n')
fprintf('                                     MÉTODO DE HEUN\n')
fprintf('                                     --------------\n')
fprintf('\n');
syms x y
d=input(' - Introduzca la ecuación diferencial        : ');
n=input(' - Introduzca la condición y(a)=b            : ');
f1=input(' - Introduzca la función de trabajo          : ');
ya=input(' - Introduzca la condición inicial           : ');
a=input(' - Introduzca el valor de a                  : ');
b=input(' - Introduzca el valor de b                  : ');
h=input(' - Introduzca el tamaño de paso h            : ');

fprintf('\n\n');

fprintf(' - La solución de la ecuación diferencial es : \n\n\n');

m = dsolve(d,n,'x');
pretty(m);

fprintf('\n\n\n');

%Condiciones para el funcionamiento de los lazos FOR
f=f1;
w(1)=ya;
i=0;
t(1)=a;
v(1)=a;
d=0;
c=0;
g=0;

%Este for obtiene y guarda todos los valores de t
%También se utiliza para evaluar la ecuación diferencial
for p=a:h:b
   d=1+d;
   t(d)=p;
   v(d)=subs(m,p);
end

%Este for se usa para contabilizar las iteraciones
for s=c:1:(d-1)
   g=1+g;
   k(g)=(g-1);
end

k3=k(end);
 
%Este for obtiene los valores aproximados de solución
fprintf('-------------------------------------------------------------------------------------------------------');
fprintf('\n');
fprintf('                                      FÓRMULAS DE CADA ITERACIÓN');
fprintf('\n');
fprintf('-------------------------------------------------------------------------------------------------------');
fprintf('\n\n');
fprintf('- w0 = %1.5f ',ya);
fprintf('\n');

for j=a:h:(b-h)
   i=1+i;
   w(i+1)=w(i)+((h/4)*(subs(f,{x,y},{t(i),w(i)})))+(((3/4)*h)*(subs(f,{x,y},{(t(i)+((2/3)*h)),(w(i)+(((2/3)*h)*(subs(f,{x,y},{t(i),w(i)}))))})));
   fprintf('\n');
   fprintf('- w%1.0f = w%1.0f + h/4 f(t%1.0f,w%1.0f) + 3/4 h f(t%1.0f + 2/3 h,w%1.0f + 2/3 h f(t%1.0f,w%1.0f))',i,i-1,i-1,i-1,i-1,i-1,i-1,i-1);
   fprintf('\n');
   fprintf('- w%1.0f = w%1.0f + %1.5f f(%1.9f,w%1.0f) + %1.5f f(%1.9f + %1.5f,w%1.0f + %1.5f f(%1.9f,w%1.0f))',i,i-1,h/4,t(i),i-1,(3/4)*h,t(i),(2/3)*h,i-1,(2/3)*h,t(i),i-1);
   fprintf('\n');
end

fprintf('\n');
fprintf('-------------------------------------------------------------------------------------------------------');
fprintf('\n');

%Presentación de los datos
fprintf('\n\n');

fprintf('              i             ti              wi+1                  y(t)');  
fprintf('\n\n');

for k1=0:k3
   k2=k1+1;
   fprintf('\n');
   fprintf('              %1.0f        %10.9f        %10.9f         %10.9f',k(k2),t(k2),w(k2),v(k2));
   fprintf('\n');                                    
end

fprintf('\n');  

 

METODO EULER MODIFICADO

% - Introduzca la ecuación diferencial        : 'Dy=y-(x^2)+1'
% - Introduzca la condición y(a)=b            : 'y(0)=0.5'
% - Introduzca la función de trabajo          : y-(x^2)+1
% - Introduzca la condición inicial           : 0.5
% - Introduzca el valor de a                  : 0
% - Introduzca el valor de b                  : 1
% - Introduzca el tamaño de paso h            : 0.1

clear all
clc
disp('---------------------------------------------------------------------------');
disp('Método modificado de euler');
disp('---------------------------------------------------------------------------');
fprintf('\n')

syms x y
d=input(' - Ingrese la ecuación diferencial    : ');
n=input(' - Ingrese la condición y(a)=b        : ');
f1=input(' - Ingrese la función de trabajo      : ');
ya=input(' - Ingrese la condicion inicial       : ');
a=input(' - Ingrese el valor de a            : ');
b=input(' - Ingrese el valor de b            : ');
h=input(' - Ingrese el valor del paso h        : ');
fprintf('\n');

fprintf(' - La solución de la ecuación diferencial es : \n\n\n');
m = dsolve(d,n,'x');
pretty(m);
fprintf('\n');  
%Condiciones para el funcionamiento de los lazos FOR
f=f1;
w(1)=ya;
i=0;
t(1)=a;
v(1)=a;
d=0;
c=0;
g=0;
%Este for obtiene y guarda todos los valores de t
%También se utiliza para evaluar la ecuación diferencial
for p=a:h:b
   d=1+d;
   t(d)=p;
   v(d)=subs(m,p);
end
%Este for se usa para contabilizar las iteraciones
for s=c:1:(d-1)
   g=1+g;
   k(g)=(g-1);
end
k3=k(end);
%Este for obtiene los valores aproximados de solución

disp('---------------------------------------------------------------------------');
disp('Formulas de cada iteracion');
disp('---------------------------------------------------------------------------');
fprintf('\n')
fprintf('- w0 = %1.5f ',ya);
fprintf('\n');
for j=a:h:(b-h)
   i=1+i;
   w(i+1)=w(i)+((h/2)*(subs(f,{x,y},{t(i),w(i)})))+((h/2)*(subs(f,{x,y},{(t(i)+h),(w(i)+(h*(subs(f,{x,y},{t(i),w(i)}))))})));
 fprintf('\n');
 fprintf('- w%1.0f = w%1.0f + h/2 f(t%1.0f,w%1.0f) + h/2 f(t%1.0f + h,w%1.0f + h f(t%1.0f,w%1.0f))',i,i-1,i-1,i-1,i-1,i-1,i-1,i-1);
  fprintf('\n');
fprintf('- w%1.0f = w%1.0f + %1.5f f(%1.9f,w%1.0f) + %1.5f f(%1.9f + %1.5f,w%1.0f + %1.5f f(%1.9f,w%1.0f))',i,i-1,h/2,t(i),i-1,h/2,t(i),h,i-1,h,t(i),i-1);
fprintf('\n');
end
fprintf('\n');
fprintf('--------------------------------------------------------------------------');
fprintf('\n');
%Presentación de los datos
fprintf('\n');
fprintf('              i             ti               wi+1                 y(t)');  
fprintf('\n');
for k1=0:k3
   k2=k1+1;
   fprintf('\n');
   fprintf('              %1.0f        %10.9f        %10.9f         %10.9f',k(k2),t(k2),w(k2),v(k2));
   fprintf('\n');                                    
end
fprintf('\n');  

METODO DEL PUNTO MEDIO

%MÉTODO DEL PUNTO MEDIO


% - Introduzca la ecuación diferencial        : 'Dy=y-(x^2)+1'
% - Introduzca la condición y(a)=b            : 'y(0)=0.5'
% - Introduzca la función de trabajo          : y-(x^2)+1
% - Introduzca la condición inicial           : 0.5
% - Introduzca el valor de a                  : 0
% - Introduzca el valor de b                  : 1
% - Introduzca el tamaño de paso h            : 0.1

fprintf('\n');
clear all
clc
fprintf('                                  ----------------------\n')
fprintf('                                  MÉTODO DEL PUNTO MEDIO\n')
fprintf('                                  ----------------------\n')
fprintf('\n');
syms x y
d=input(' - Introduzca la ecuación diferencial        : ');
n=input(' - Introduzca la condición y(a)=b            : ');
f1=input(' - Introduzca la función de trabajo          : ');
ya=input(' - Introduzca la condición inicial           : ');
a=input(' - Introduzca el valor de a                  : ');
b=input(' - Introduzca el valor de b                  : ');
h=input(' - Introduzca el tamaño de paso h            : ');

fprintf('\n\n');

fprintf(' - La solución de la ecuación diferencial es : \n\n\n');

m = dsolve(d,n,'x');
pretty(m);

fprintf('\n\n\n');

%Condiciones para el funcionamiento de los lazos FOR
f=f1;
w(1)=ya;
i=0;
t(1)=a;
v(1)=a;
d=0;
c=0;
g=0;

%Este for obtiene y guarda todos los valores de t
%También se utiliza para evaluar la ecuación diferencial
for p=a:h:b
   d=1+d;
   t(d)=p;
   v(d)=subs(m,p);
end

%Este for se usa para contabilizar las iteraciones
for s=c:1:(d-1)
   g=1+g;
   k(g)=(g-1);
end

k3=k(end);
 
%Este for obtiene los valores aproximados de solución
fprintf('      ----------------------------------------------------------------------------');
fprintf('\n');
fprintf('                               FÓRMULAS DE CADA ITERACIÓN');
fprintf('\n');
fprintf('      ----------------------------------------------------------------------------');
fprintf('\n\n');
fprintf('       - w0 = %1.5f ',ya);
fprintf('\n');

for j=a:h:(b-h)
   i=1+i;
   w(i+1)=w(i)+(h*(subs(f,{x,y},{(t(i)+h/2),w(i)+((h/2)*(((subs(f,{x,y},{t(i),w(i)})))))})));  
   fprintf('\n');
   fprintf('       - w%1.0f = w%1.0f + h f(t%1.0f + h/2,w%1.0f + h/2 f(t%1.0f,w%1.0f))',i,i-1,i-1,i-1,i-1,i-1);
   fprintf('\n');
   fprintf('       - w%1.0f = w%1.0f + %1.5f f(%1.9f + %1.5f,w%1.0f + %1.5f f(%1.9f,w%1.0f))',i,i-1,h,t(i),(h/2),i-1,(h/2),t(i),i-1);
   fprintf('\n');
end

fprintf('\n');
fprintf('      ----------------------------------------------------------------------------');
fprintf('\n');

%Presentación de los datos

fprintf('\n\n');

fprintf('              i             ti               wi+1                 y(t)');  

fprintf('\n\n');

for k1=0:k3
   k2=k1+1;
   fprintf('\n');
   fprintf('              %1.0f        %10.9f        %10.9f         %10.9f',k(k2),t(k2),w(k2),v(k2));
   fprintf('\n');                                    
end

fprintf('\n');  

METODO DE EULER

%MÉTODO DE EULER

% - Introduzca la ecuación diferencial        : 'Dy=cos(2*x)+sin(3*x)'
% - Introduzca la condición y(a)=b            : 'y(0)=1'
% - Introduzca la función de trabajo          : cos(2*x)+sin(3*x)
% - Introduzca la condición inicial           : 1
% - Introduzca el valor de a                  : 0
% - Introduzca el valor de b                  : 1
% - Introduzca el tamaño de paso h            : 1/4

fprintf('\n');
clear all
clc
fprintf('                                     ---------------\n')
fprintf('                                     MÉTODO DE EULER\n')
fprintf('                                     ---------------\n')
fprintf('\n');
syms x y
d=input(' - Introduzca la ecuación diferencial        : ');
n=input(' - Introduzca la condición y(a)=b            : ');
f1=input(' - Introduzca la función de trabajo          : ');
ya=input(' - Introduzca la condición inicial           : ');
fprintf('Introduzca el intervalo de evaluacion.\n');
a=input(' - Desde                                     : ');
b=input(' - Hasta                                     : ');
h=input(' - Introduzca el tamaño de paso h            : ');

fprintf('\n\n');

fprintf(' - La solución de la ecuación diferencial es : \n\n\n');

m = dsolve(d,n,'x');
pretty(m)

fprintf('\n\n\n');

%Condiciones para el funcionamiento de los lazos FOR
f=f1;
w(1)=ya;
i=0;
t(1)=a;
v(1)=a;
d=0;
c=0;
g=0;

%Este for obtiene y guarda todos los valores de t
%También se utiliza para evaluar la ecuación diferencial
for p=a:h:b
   d=1+d;
   t(d)=p;
   v(d)=subs(m,p);
end

%Este for se usa para contabilizar las iteraciones
for s=c:1:(d-1)
   g=1+g;
   k(g)=(g-1);
end

k3=k(end);

%Este for obtiene los valores aproximados de solución

%También imprime en pantalla la fórmula de la ecuación para
%cada iteración

fprintf('           ------------------------------------------------');
fprintf('\n');
fprintf('                      FÓRMULAS DE CADA ITERACIÓN');
fprintf('\n');
fprintf('           ------------------------------------------------');
fprintf('\n\n');
fprintf('            - w0 = %10.15f ',ya);
fprintf('\n');

for j=a:h:(b-h)
   i=1+i;
   w(i+1)=w(i)+(h*(subs(f,{x,y},{t(i),w(i)})));  
   fprintf('\n');
   fprintf('            - w%1.0f = w%1.0f + h f(t%1.0f,w%1.0f)',i,i-1,i-1,i-1);
   fprintf('\n');
   fprintf('            - w%1.0f = w%1.0f + %1.5f f(%10.15f,%10.15f)',i,i-1,h,t(i),w(i));
   fprintf('\n');
end

fprintf('\n');
fprintf('           ------------------------------------------------');
fprintf('\n');

fprintf('\n\n');

fprintf('      i               ti                 wi+1                   y(t)');  

fprintf('\n\n');

for k1=0:k3
   k2=k1+1;
   fprintf('\n');
   fprintf('      %1.0f        %10.15f        %10.15f         %10.15f',k(k2),t(k2),w(k2),v(k2));
   fprintf('\n');                                    
end

fprintf('\n');




 

Ecuaciones Diferenciales de primer orden

·         Archivo M

function  dy=g1(t,y)

dy=funcion en términos t y Y;

·         Para cualquier tipo de ode:

[t y]=ode113(@g1,[valor inferior de t:espaciado h:valor inferior de t],valor inicial];

%En los cmondos anteriores puede ser cualquier tipo de ode

·         Para la solución exacta

Exacta=dsolve(‘Dy=funcion en términos de y y t’,’Y()=valor inicial’)

Y=subs(Exacta,t)

·         Error

Error=norm(Y-y,inf)

Ecuación de orden superior de 2 ecuaciones  

%En este caso se hará un proceso a través de método del anulador  para obtener la funcion DY Y DX

Archivo M

function p=f1(t,w)

p=zeros(size(w));

p=[numerox 0;0 numeroy]*w+[funcion restante de dx;function restante de dy];

·         Para cualquier tipo de Ode

t=[valor inferior:espaciado h;valor superior];

[t xy]=ode45(@f1,t,[valor X inicial;valor Y inicial])

·         Exacto

[x y]=dsolve(‘funcion DX','Funcion DY','x(0)=0','y(0)’)

ExactoX=subs(x,t);

ExactoY=subs(y,t);

Exacto=[ExactoX ExactoY]

errorx=norm(ExactoX-xy(:,1));

errory=norm(ExactoY-xy(:,2));

error=[errorx errory]

Comados para Integrales

Resolucion de la primera intregral

·         Para quad, quadl, quadv, quadgk:

Los comandos se explicaran con base a un ejemplo:

Archivo M

function  f=g1(x) %Donde g1 puede ser cualquier nombre.

f=exp(x).*cos(5*x); %Se ingresa cualquier tipo de funcion.

%Cuando se ingresan las funciones en el archivo M se debe de pones un punto cuando la variable x se divide o se multiplica con otra función o variable x.

 quad

PrimeraIntegral=quad(@g1,limite inferior (a), limite superior (b))

quadl

PrimeraIntegral=quadl(@g1,limite inferior (a), limite superior (b))

quadv

PrimeraIntegral=quadv(@g1,limite inferior (a), limite superior (b))

quadgk

PrimeraIntegral=quadgk(@g1,limite inferior (a), limite superior (b))

%Donde @g1 se puede cambiar según el nombre dado al Archivo M.

%Exacto

Exacto=double(int(‘funciona integrar sin los puntos de separación’,’x’,limite inferior(a),limite superior(b)))

·         Para trapz

x=linspace(limite inferior(a),limite superior(b))

y=función en términos de x y separados por puntos en división y multiplicación;

Integral=trapz(x,y)

  

Resolución de la segunda integral

·         Uso del dblquad:

Archivo M

function f=g7(x,y)

f=Fucion en terminos de ‘x’ y ‘y’;

dobleintegral=dblquad(@g7,intervalo inferio de (dx), intervalo superior de (dx), intervalo inferio de (dy), intervalo superior de (dy));

·         Uso de int

Ejemplo: si el orde de integración es la siguiente dxdy y ax=intervalo inferior de x, bx=intervalo superior de x, ay=intervalo inferior de y, by=intervalo superior de y.

integralconINT=int(int(‘funcion en términos X y Y’,’x’,ax,bx),’y’,ay,by)

Resolución de la tercera  integral

·         Uso del triplequad:

Archivo M

function f=g9(x,y,z)

f=Fucion en terminos de ‘x’ y ‘y’;

tripleintegral=triplequad(@g9,intervalo inferio de (dx), intervalo superior de (dx), intervalo inferio de (dy), intervalo superior de (dy), intervalo inferio de (dz), intervalo superior de (dz));

·         Uso de int

Ejemplo: si el orde de integración es la siguiente dxdydz y ax=intervalo inferior de x, bx=intervalo superior de x, ay=intervalo inferior de y, by=intervalo superior de y, az=intervalo inferior de z, bz=intervalo superior de z.

integralconINT=int(int(int(‘funcion en términos X , Yy Z’,’x’,ax,bx),’y’,ay,by),’z’,az,bz)

Comandos para derivadas.

Derivadas

·         Derivada solo con datos con “h” uniforme:

%Ingreso de los datos de X:

X=[ingreso de datos con espacios];

%Ingreso de los datos de Y:

Fx=[ingreso de datos con espacios];

PrimeraDerivada=diff(fx)./h

%Donde h es el valor de separación de cada dato de X

·         Derivada solo con datos con “h” NO uniforme:

%Ingreso de los datos de X:

X=[ingreso de datos con espacios];

%Ingreso de los datos de Y:

Fx=[ingreso de datos con espacios];

PrimeraDerivada=diff(fx)./diff(fx)

·         Primera derivada con función

syms x

P=[datos de x con espacios];

Y=[datos de fx con espacios];

F=función en términos de x;

%Derivada aproximada

Aprox=diff(y)./diff(P)

%Derivada exacta

Df=diff(f);    %se va a derivar la función

Exact=subs(Df,x,P(1:end-1))

%Error

error=norm(Exact-Aprox,inf)

·         Segunda derivada sin datos:

X=[Ingreso de datos de x con espacios];

Y=[ingreso de los valores fx con espacios];

Primeraderivada=diff(Y)./diff(X);

SegundaDerivada=diff(Primeraderivada)./diff(X(1:end-1))

·         Segunda derivada con función:

syms x

P=[Ingreso de datos de x con espacios];

Y=[ingreso de los valores fx con espacios];

F=ingreso de la función en términos de x;

PrimeraDerivadaAprox=diff(y)./diff(P);

SegundaDerivadaAprox=diff(PrimeraDerivadaExacta)./diff(P(1:end-1))

%Exacta

Df=diff(f,x); %Primera Derivada

Dff=diff(Df,x); %Segunda Derivada

Segundaderivadaexacta=subs(Dff,x,P(1:end-2))

%Error

Error=nrom(Segundaderivadaexacta- PrimeraDerivadaAprox,inf)

Formulas y factores de conversión

• Temperatura.          

      De F a K -------> (5/9)*(F - 32)+273.15

         De F a C------->(5/9)*(F - 32)

         De C a F-------->(9/5)*C + 32

•  Masa.

         1kg = 2.2 lb

     

•     Unidades Metricas

         1m = 100cm

         1Km = 1000m

         1pul = 1/12pies

           

Interpolacion y Extrapolacion a Travez de "interp1(..)"

• Cubica sin función

Interpolacion:

X1=[Ingreso de datos de x];

Y=[Ingreso de los espejos Y];

interp1(X1,Y,valor a interpolar,'cubic')

Extrapolacion:

X1=[Ingreso de datos de x];

Y=[Ingreso de los espejos Y];

interp1(X1,Y,valor a interpolar,'cubic','extrap')

• Cubica con función 

Interpolacion: 

syms x

X1=[Ingreso de datos de x];

f=función para interpolar en terminos de "x";

Y=subs(f,X1);

interp1(X1,Y,valor a interpolar,'cubic')

Extrapolacion:

syms x

X1=[Ingreso de datos de x];

f=función para interpolar en terminos de "x";

Y=subs(f,X1);

interp1(X1,Y,valor a interpolar,'cubic','extrap')

Interpolacion y Extrapolacion a Travez de "interp1(..)"

• Lineal sin función

Interpolacion:

X1=[Ingreso de datos de x];

Y=[Ingreso de los espejos Y];

interp1(X1,Y,valor a interpolar,'linear')

Extrapolacion:

X1=[Ingreso de datos de x];

Y=[Ingreso de los espejos Y];

interp1(X1,Y,valor a interpolar,'linear','extrap')

• Lineal con función 

Interpolacion: 

syms x

X1=[Ingreso de datos de x];

f=función para interpolar en terminos de "x";

Y=subs(f,X1);

interp1(X1,Y,valor a interpolar,'linear')

Extrapolacion:

syms x

X1=[Ingreso de datos de x];

f=función para interpolar en terminos de "x";

Y=subs(f,X1);

interp1(X1,Y,valor a interpolar,'linear','extrap')

Interpolacion y Extrapolacion a Travez de "interp1(..)"

• Escalonada sin función

Interpolacion:

X1=[Ingreso de datos de x];

Y=[Ingreso de los espejos Y];

interp1(X,Y,valor a interpolar,'spline')

Extrapolacion:

X1=[Ingreso de datos de x];

Y=[Ingreso de los espejos Y];

interp1(X,Y,valor a interpolar,'spline','extrap')

• Escalonada con función 

Interpolacion: 

syms x

X1=[Ingreso de datos de x];

f=función para interpolar en terminos de "x";

Y=subs(f,X1);

interp1(X,Y,valor a interpolar,'spline')

Extrapolacion:

syms x

X1=[Ingreso de datos de x];

f=función para interpolar en terminos de "x";

Y=subs(f,X1);

interp1(X,Y,valor a interpolar,'spline','extrap')

Interpolacion de Diferencias Divididas con Datos y Función 

clear all
clc
syms x
disp('MÉTODO DE LAS DIFERENCIAS DIVIDIDAS')
fprintf('-------------------------------------\n');
val=input('Valor a interpolar x: ');
dato=input('Datos Ordenados[X0 X1 X2 ... Xn]: ');
t=length(dato);
fprintf('Valores de la función:\n\t1-Utilizar una función.\n\t2-Ingresar valores\n');
opc=input('Opción: ');
switch opc
    case 1
        f=input('Función f(x): ');
        fun=zeros(1,t);
        fprintf('Valores de F(x): ');
        for i=1:t
            fun(i)=subs(f,dato(i));
            fprintf('%8.09f, ', fun(i));
        end
    case 2
        fun=input('Valores F(x) [F(X0) F(X1) ... F(Xn)]: ');
end
DD=zeros(t);
fprintf('\nObteniendo las diferencias divididas:');
DD(:,1)=fun';
%for i=1:t
%    DD(i,1)=fun(i);
%end
for j=2:t
    fprintf('\n - Columna %2.0f de diferencias divididas\n',j);
    for i=j:t
        DD(i,j)=(DD(i,j-1)-DD(i-1,j-1))/(dato(i)-dato(i-j+1));
        fprintf('(%8.15f-(%8.15f))/(%3.3f-(%3.3f)) = %8.09f\n',DD(i,j-1),DD(i-1,j-1),dato(i),dato(i-j+1),DD(i,j));
    end
end
DD %#ok<NOPTS>
fprintf('\nPolinomio grado %1.0f\n',t-1);
pol=DD(1,1);
fprintf('P%1.0f(x)= A0 + ',t-1);
for i=2:t
    fprintf('A%1.0f',i-1);
    for j=1:i
        fprintf('(x-x%1.0f)',j-1);
    end
    if i==t
        fprintf('\n');
    else
        fprintf(' + ');
    end
end
fprintf('P%1.0f(x)= %8.09f + ',t-1,DD(1,1));
for i=2:t
    factor=DD(i,i);
    fprintf('(%8.09f)*',DD(i,i));
    for j=1:i-1
        fprintf('(x-(%3.3f))',dato(j));
        factor=factor*(val-dato(j));
    end
    pol=pol+factor;
    if i==t
        fprintf('\n');
    else
        fprintf(' + ');
    end
end
fprintf('P%1.0f(%3.3f)= %9.15f\n\n',t-1,val,pol);
if opc==1
    fprintf('Valor Exacto de la Función: %9.15f',subs(f,val));
    fprintf('\nError: %e\n\n',abs(pol-subs(f,val)));

end

Interpolacion de Neville con Datos Y Función 

syms x;
vx = input('Ingrese el vector [X0 X1 X2 ......Xn]: ');
xint = input('Ingrese el punto a interpolar: ');
disp('Seleccione como trabajara el metodo');
opcion = input('1) Vector;  2) Funcion;  opcion: ');
vy=zeros(length(vx));
if(opcion == 1)
    vy = input('Ingrese el vector [f(x0) f(x1) f(x2).....f(xn)]: ');
elseif(opcion == 2)
    f = input('Ingrese la funcion de trabajo: ');
    for i = 1:length(vx)
        vy(i) = subs(f, vx(i));
    end
else
    disp('Opcion invalida :( .|.\n');
end
Q = zeros(length(vx), length(vx));
for i = 1 : length(vx)
    Q(i, 1) = vy(i);
    fprintf('Q(%1.0f,0) = %3.15f\n', i-1, vy(i));
end
disp('-------');
for j = 2:length(vx) %x
    for i = j: length(vx) %y
        Q(i, j) = ( (xint - vx(i-j+1))*Q(i, j-1) - (xint - vx(i) )*Q( i-1, j-1) ) / (vx(i) - vx(i-j+1));
        fprintf('Q(%1.0f,%1.0f) = ((x-x%1.0f)*Q(%1.0f,%1.0f)-(x-x%1.0f)*Q(%1.0f,%1.0f))/(x%1.0f-x%1.0f) = %3.15f \n',i-1,j-1, i-j, i-1, j-2, i-1, i-2, j-2, i-1, i-j, Q(i,j));%a todos se les rest -1 y lo demaa squeda igyal
    end
    disp('-------');
end
fprintf('La respuesta es: %3.15f \n', Q(length(vx), length(vx)));

Interpolacion de Hermite con Diferencias Divididas con Datos y Función 

clc;
syms x
disp('INTERPOLACION DE HERMITE POR DIFERENCIAS DIVIDIDAS')
fprintf('----------------------------------------------------------\n');
val=input('Valor a interpolar x: ');
vars=input('Datos Ordenados[X0 X1 X2 ... Xn]: ');
t=2*length(vars);
dato=zeros(1,t);
for i=1:t/2
    dato(2*i-1)=vars(i);
    dato(2*i)=vars(i);
end
fprintf('Valores de la función:\n\t1-Utilizar una función.\n\t2-Ingresar valores\n');
opc=input('Opción: ');
switch opc
    case 1
        f=input('Función f(x): ');
        fundif=zeros(1,t/2);
        fun=zeros(1,t/2);
        fprintf('Valores de F(x): ');
        for i=1:t/2
            fun(i)=subs(f,vars(i));
            fprintf('%0.15f,', fun(i));
        end
        fprintf('\nValores de F´(x): ');
        for i=1:t/2
            fundif(i)=subs(diff(f,x),vars(i));
            fprintf('%f,', fundif(i));
        end
    case 2
        fun=input('Valores F(x) [F(X0) F(X1) ... F(Xn)]: ');
        fundif=input('Valores F´(x) [F´(X0) F´(X1) ... F´(Xn)]: ');
end
DD=zeros(t);
fprintf('\nObteniendo las diferencias divididas:');
for i=1:t/2
    DD(2*i-1,1)=fun(i);
    DD(2*i,1)=fun(i);
end
for j=2:t
    for i=j:t
        if dato(i)~=dato(i-j+1)
            DD(i,j)=(DD(i,j-1)-DD(i-1,j-1))/(dato(i)-dato(i-j+1));   
        else
            DD(i,j)=fundif(i/2);
        end
    end
end
DD %#ok<NOPTS>
fprintf('\nPolinomio grado %1.0f\n',t-1);
h=DD(1,1);
fprintf('H%1.0f(x)= A0 + ',t-1);
for i=2:t
    fprintf('A%1.0f',i-1);
    for j=1:i
        fprintf('(x-x%1.0f)',j-1);
    end
    if i==t
        fprintf('\n');
    else
        fprintf(' + ');
    end
end
fprintf('H%1.0f(x)= %8.15f + ',t-1,DD(1,1));
for i=2:t
    factor=DD(i,i);
    fprintf('%8.15f*',DD(i,i));
    for j=1:i-1
        fprintf('(x-(%3.15f))',dato(j));
        factor=factor*(val-dato(j));
    end
    h=h+factor;
    if i==t
        fprintf('\n');
    else
        fprintf(' + ');
    end
end
fprintf('H%1.0f(%3.15f)= %9.15f\n\n',t-1,val,h);
if opc==1
    fprintf('Valor Exacto de la Función: %9.15f',subs(f,val));
    fprintf('\nError: %e\n\n',abs(h-subs(f,val)));

end