[Matlab] Metodo a due passi per Cauchy

[5mrkv]

L'equazione che ho usato è
\[
\begin{cases}
y'(t)&=4ty^{1/2}(t) \\
y(0)&=1 \\
\end{cases}
\]
Con soluzione
\[
y(t)=(t^{2}+1)^{2}
\]
Il metodo da utilizzare per risolverla è
\[
\begin{cases}
a&=f(t_{i+1},y_{i+1}) \\
b&=f(t_{i},y_{i}) \\
y_{n+2}&=(1+\alpha)f_{n+1}-\alpha f_{i}+(h/2)[(3-\alpha)a-(1+\alpha)b] \\
\end{cases}
\]
L'esercizio chiede di scrivere una funzione che implementi il metodo verificandone il funzionamento per \(\alpha=0,-5\). Il problema è che se con \(\alpha=0\) funziona, non funziona con \(\alpha=-5\). Nel programma ho usato prima eulero per calcolare l'approssimazione al primo passo.

f.m

function a=f(t,y)
         a=4*t*y^(1./2);

multistep.m

function [t ris]=multistep(f,ti,tf,h,yi,k)

%inizializzo
t=ti:h:tf;
y_app(1)=yi;

c1=(1+k);
c2=(3-k);
c3=(1+k);

%metodo eulero
x1=f(t(1),y_app(1));
y_app(2)=y_app(1)+h*x1;

%metodo multistep
for i=1:length(t)-2,
    x2=f(t(i+1),y_app(i+1));
    x3=f(t(i),y_app(i));
    x4=(c2*x2-c3*x3);
    y_app(i+2)=c1*y_app(i+1)-k*y_app(i)+(h./2)*x4;
end

ris=y_app;
sol=power((power(t,2)+1),2);

%confronto i grafici
plot(t,ris,'r-.',t,sol,'g-');
legend('Approssimata','Esatta');
title('Metodo multistep a due passi');
xlabel('Tempo');
ylabel('y*(t), y(t)');

script.m

ti=0;
tf=2;
h=0.1;
yi=1;
k=0;

[t,ris]=multistep(@f,ti,tf,h,yi,k);

L'errore è

Warning: Imaginary parts of complex X and/or Y arguments
ignored 
> In multistep at 27
  In script at 8

Si vede che vengono fuori delle radici complesse. Avete qualche idea? Qualcuno conosce questo metodo?

[5mrkv]

bamp