syms x;
x0=0;       %% Anfangswert y(x0)=y0
y0=1;
p=x*tan(x);            %% p-Funktion in x 
q=(x-1)*sin(x);        %% q-Funktion der DGL y' + p * y = q

n=7;        %% Grad des Ergebnispolynoms

p=p+x-x;
q=q+x-x;
a=zeros(1,n+1);
b=zeros(1,n+1);
c=zeros(1,n+1);
a(1)=y0;
b(1)=subs(p,x,x0);
c(1)=subs(q,x,x0);

for i=2:n+1
    b(i)=subs(diff(p,i-1)/factorial(i-1),x,x0);
    c(i)=subs(diff(q,i-1)/factorial(i-1),x,x0);
end


for i=2:n+1
    res=0;
    for j=1:i-1
        res=res + a(i-j)*b(j);
    end
    a(i)=(c(i-1)-res)/(i-1);
end

syms poly

poly=0;
for i=1:n+1
    poly = poly + a(i)*(x-x0)^(i-1); 
end

vpa(poly,5)

val=[];
err=[];
for xi=x0:0.01:x0+1 
    val=[val,subs(poly,x,xi)];
    err=[err,subs(diff(poly,x)+p*poly-q,x,xi)];
end

figure(1)
plot(x0:0.01:x0+1,val); %% Ergebnis zwischen x0 und x0+1 plotten