function y =desapotenciasn(a0,a1,N)
%C.S. y N. en ED, curso 0708
%function y =desapotenciasn(a0,a1,n)
%resuelve un problema de Cauchy para la ecuacion 
%%y''+exp(x)*y'+(1+x^2)y=0 con condiciones inciales y(0)=a0, y'(0)=a1
%%calculando los N primeros terminos del desarrollo en serie de potencias
%de la solucion. Lee los datos iniciales a0, a1 y N
%inicializamos un vector
a=zeros(N,1);
c=zeros(N,1);
factor=zeros(N,1);
a(1)=a1;
a(2)=-(a(1)+a0)/2;
c(1)=a(1)+2*a(2);
a(3)=-(a(1)+c(1))/6;
c(2)=3*a(3)+2*a(2)+a(1)/2;
a(4)=-(a0+a(2)+c(2))/12;
for i=3:N-2
    z=0;
    for j=1:i
        z=z+a(i+1-j)*(i-j+1)/factorial(j);
    end
    c(i)=z+a(i+1)*(i+1);
    a(i+2)=-((a(i-2)+a(i)+c(i))/((i+2)*(i+1)));
end
syms x;
h=a0;
for i=1:N
    h=h+a(i)*x^i;
end
nterminosdesarrollo=h
ezplot(h,0,4)
