% grxxm063.m
% Taylorreihen-Methode fuer  y' = - 2 t y^2, t > 0, y(0) = y0 > 0
% Loesung des AWP:  y(t) = y0/(y0 t^2 + 1) 
% y_{n+1} = y_n + c_1 h + c_2 h^2 + c_3 h^3 + c_4 h^4

% Anfangswert
y0=1;   %%% Eingabewert %%%

% Integrationsintervall [0,T]
T=1;

% Schrittweite h
n=20;   %%% Eingabewert %%%
h=T/n;
h2=h^2;
h3=h^3;
h4=h^4;

T=zeros(n+1,7);

% Belegung von T(1,1) und T(1,2)
T(1,1)=0;
T(1,2)=y0;


% Iteration
for k=1:n
    T(k+1,1)=k*h;
    T(k,3)=-2*T(k,1)*T(k,2)^2;
    T(k,4)=-(T(k,2)+2*T(k,3)*T(k,1))*T(k,2);
    T(k,5)=-(4*T(k,3)*T(k,2)+2*T(k,1)*(T(k,3)^2+2*T(k,4)*T(k,2)))/3;
    T(k,6)=-(T(k,3)^2/2+T(k,4)*T(k,2)+T(k,1)*(T(k,3)*T(k,4)+T(k,5)*T(k,2)));
    T(k+1,2)=T(k,2)+T(k,3)*h+T(k,4)*h2+T(k,5)*h3+T(k,6)*h4;
    T(k+1,7)=T(k+1,2)-1/(T(k+1)^2+1/y0);
end;

% Tabelle 
disp('-------------------------------------------------------------------------');
disp('           Taylorreihen-Methode fuer  dy/dt = - 2 t y^2');
disp('-------------------------------------------------------------------------');
disp('t_n      y_n       c_1(n)     c_2(n)     c_3(n)     c_4(n)     Fehler ');
disp('-------------------------------------------------------------------------');

for s=1:n+1
    fprintf('%2.3f %9.6f %10.6f %10.6f %10.6f %10.6f %13.10f',...
        T(s,1),T(s,2),T(s,3),T(s,4),T(s,5),T(s,6),T(s,7));
    fprintf('\n');
end;
    




