1. Algoritmo para el método de integración de Euler para la ecuación diferencial 
y' =  x y. Comparación con el método numérico NDSolve de Mathematica

Clear["@"];
f[x_, y_] :=  x y;
(*Consideramos el problema de valores iniciales   
 y' = f (x, y),  y (a) = (Subscript[y, 1]) = :ci,   (VI  .5  .1)
y queremos encontrar una solución aproximada en el  intervalo [a, b] .  Para ello comenzamos dividiendo el
intervalo [a, b] en N = 10 ene subintervalos iguales de longitud (paso) h = (b - a)/N,*)
a = 
        0;
b = 
        0.5;
ci = 1;
h = (b - a)/(10 ene);
  
x[1] = a;
x[n_] := x[n - 1] + h;
(* Para valores de h pequeños, \(la\ función\ y \((x)\)\ en\ el\ intervalo\ [a, 
                x\_2]\ se\ puede\ aproximar\ bien\ por\ la

            
            función\ cuya\ gráfica\ es\ la\ recta\ tangente\ en\ a\), y (x) = y (a) + y' (a) (x - a) = 
                y (a) + f (a, (Subscript[y, 1])) (x - a), \(donde\ hemos\ sustituido\ y' \((a)\)\ por\ la\ expresión\ que\ tiene\ en
\(nuestro\ problema\ de\ valores\ iniciales . 
                  Por\ lo\ tanto\)\),aproximadamente,\(y\_2 = 
              y \((x\_2)\)\ ~\ y \((a)\) + 
                f \((a, y\_1)\)\ h . 

                    A\ partir\ de\ este\ valor\ aproximado\ de\ y \((x\_2)\)\),\(podemos

            
            usar\ el\ mismo\ argumento\ para\ calcular\ el\ valor\ aproximado

            
            de\ y \((x\_3)\) : \ 
              y \((x\_3)\)\ ~\ y \((x\_2)\)\  + \ f \((x\_2, y\_2)\)\ h\), y así sucesivamente .  
                  Escribiremos todos los pasos a hacer así :   (Subscript[y, n]) := y ((Subscript[x, n])) ~ y ((Subscript[x, n - 1]))  + 
                f ((Subscript[x, n]), (Subscript[y, n])) h,   1 ≤ n ≤ N.    (VI  .5  .2)
Esta es la fórmula del método de Euler o método de la
tangente,que aproxima la solución de (VI  .5  .1) por (VI  .5  .2)*)
y[1] = 
        ci;
y[n_] := y[n - 1] + f[x[n - 1], y[n - 1]] h;
  
Eulerfacire = 
        Table[ListPlot[Table[{x[n], y[n]}, {n, 1, 10 ene + 1}], 
            PlotJoined -> True,  
            PlotStyle -> {Hue[0.3 ene]}], {ene, 1, 
            2}];
  Eulerfacipu = 
          Table[ListPlot[Table[{x[n], y[n]}, {n, 1, 10 ene + 1}],  
              PlotStyle -> {Hue[0.3 ene]}], {ene, 1, 
              2}];
  ndsolveypf = 
          Plot[Evaluate[
              z[x]  /. NDSolve[{z'[x] ==  x z[x],  
                    z[a] == ci}, z, {x, a, b}]], {x, a, b}, 
            PlotStyle -> {Hue[0.9]}];
  
Show[Eulerfacire, 
        Eulerfacipu, ndsolveypf, 
        DisplayFunction -> $DisplayFunction];

Solución  del problema de valores iniciales 
                                 y' = cos(15 x y),   y(0) = 1, 
en el intervalo [ 0, 0.5 ] usando el desarrollo de Taylor . Comparación con el método numérico NDSolve de Mathematica

y[n_, x_] := ci + (∑_(i = 1)^n) c[i] (x)^(i)

g[n_, x_] := Normal[Series[f[x, y[n, x]], {x, 0, n - 1}]]

yp[n_, x_] := D[y[n, x], x];

coef[n_] := 
      Solve[CoefficientList[yp[n, x], x] == 
          CoefficientList[g[n, x], x]]

f[x_, y_] := Cos[15 x y];

ci = 1;

a = 0;  b = 0.5;

coef[5]

coefic[n_] := Sort[coef[n][[1]]]

coefic[5] // MatrixForm

coefic[5][[3, 2]]

c[i_, n_] := coefic[n][[i, 2]]

ysol[n_, x_] := 
      ci  +  (∑_(i = 1)^n) c[i, n] (x)^(i)

sol5 = ysol[5, x]

d5 = 
        Plot[sol5, {x, 0, 0.5}, PlotStyle ->  {Hue[0.4]}];

sol7 = ysol[7, x]

d7 = 
        Plot[sol7, {x, 0, 0.5}, PlotStyle ->  {Hue[0.6]}];

sol8 = ysol[8, x]

d8 = 
        Plot[sol8, {x, 0, 0.5}, PlotStyle ->  {Hue[0.9]}];

g1 = 
Plot[
          Evaluate[
            z[x]  /. NDSolve[{z'[x] == Cos[15 x z[x]],  
                  z[a] == 1}, z, {x, a, b}]], {x, a, b}];

Show[g1, d5, d7,  d8];

1. Algoritmo para el método de integración de Euler para la ecuación diferencial 
y' = cos(15 x y). Comparación con el método numérico NDSolve de Mathematica

Clear["@"];
f[x_, y_] := Cos[15 x y];
(*Consideramos el problema de valores iniciales   
 y' = f (x, y),  y (a) = (Subscript[y, 1]) = :ci,   (VI  .5  .1)
y queremos encontrar una solución aproximada en el  intervalo [a, b] .  Para ello comenzamos dividiendo el 
intervalo [a, b] en (N  =  10 ene) subintervalos iguales de longitud (paso) h=(b - a)/N,*)
a = 
        0;
b = 
        0.5;
ci = 1;
h = (b - a)/(10 ene);
  
x[1] = a;
x[n_] := x[n - 1] + h;
(* Para valores de h pequeños, \(la\ función\ y \((x)\)\ en\ el\ intervalo\ [a, 
                x\_2]\ se\ puede\ aproximar\ bien\ por\ la

            
            función\ cuya\ gráfica\ es\ la\ recta\ tangente\ en\ a\), y (x) = y (a) + y' (a) (x - a) = 
                y (a) + f (a, (Subscript[y, 1])) (x - a), \(donde\ hemos\ sustituido\ y' \((a)\)\ por\ la\ expresión\ que\ tiene\ en
\(nuestro\ problema\ de\ valores\ iniciales . 
                  Por\ lo\ tanto\)\),aproximadamente,\(y\_2 = 
              y \((x\_2)\)\ ~\ y \((a)\) + 
                f \((a, y\_1)\)\ h . 

                    A\ partir\ de\ este\ valor\ aproximado\ de\ y \((x\_2)\)\),\(podemos

            
            usar\ el\ mismo\ argumento\ para\ calcular\ el\ valor\ aproximado

            
            de\ y \((x\_3)\) : \ 
              y \((x\_3)\)\ ~\ y \((x\_2)\)\  + \ f \((x\_2, y\_2)\)\ h\), y así sucesivamente .  
                  Escribiremos todos los pasos a hacer así :   (Subscript[y, n]) := y ((Subscript[x, n])) ~ y ((Subscript[x, n - 1]))  + 
                f ((Subscript[x, n]), (Subscript[y, n])) h,   1 ≤ n ≤ N.    (VI  .5  .2)
Esta es la fórmula del método de Euler o método de la
tangente,que aproxima la solución de (VI  .5  .1) por (VI  .5  .2)*)
y[1] = 
        ci;
y[n_] := y[n - 1] + f[x[n - 1], y[n - 1]] h;
  
Eulerfacire = 
        Table[ListPlot[Table[{x[n], y[n]}, {n, 1, 10 ene + 1}], 
            PlotJoined -> True,   PlotStyle -> {Hue[0.3 ene]}, 
            DisplayFunction -> Identity], {ene, 1, 
            2}];
  Eulerfacipu = 
          Table[ListPlot[
              Table[{x[n], y[n]}, {n, 1, 
                  10 ene + 1}],   PlotStyle -> {Hue[0.3 ene]}, 
              DisplayFunction -> Identity], {ene, 1, 
              2}];
  ndsolveypf = 
          Plot[Evaluate[
              z[x]  /. NDSolve[{z'[x] == Cos[15 x z[x]],  
                    z[a] == ci}, z, {x, a, b}]], {x, a, b},  
            PlotStyle -> {Hue[0.9]}, 
            DisplayFunction -> Identity];
  
Show[Eulerfacire, 
        Eulerfacipu, ndsolveypf, 
        DisplayFunction -> $DisplayFunction];