Problema matlab

[deborolaroccia]

Ciao, ho il seguente problema differenziale ai valori iniziali
\(\displaystyle y′1 (t) = −2y1(t)y2(t)\)
\(\displaystyle y′2 (t) = y1(t)^2 − y2(t)^2 + y3(t)^2 − 1\)
\(\displaystyle y′3 (t) = −2(y1(t) + y2(t))y3(t)\)
\(\displaystyle y1(0) = y2(0) = y3(0) = 1 \)
Con tableau

0		1
1/2		1/2

Devo calcolare numericamente l'ordine.

Io ho scritto il seguente codice

clear all
close all
clc

mrif = 1000;
m = mrif;
A = [0,0;1/2,1/2];
b = [1/2,1/2];
c = [0,1];
tspan = [0,1];
Y = zeros(3,m+1); % Matrice che memorizza le iterazioni
y0=[1;1;1];
Y(:,1) = y0;

odefun = @(x) [-2*x(1)*x(2);(x(1))^2-(x(2))^2+(x(3))^2 -1;-2*(x(1)+x(2))*x(3)];
odeJf = @(x) [-2*x(2),-2*x(1),0; 2*x(1),-2*x(2),2*x(3);-2*x(3),-2*x(3),-2*(x(1)+x(2))];
k = 2^(-12);
m = floor((tspan(2)-tspan(1))/k)+1;
tout = zeros(1,m+1);
tout(1) = tspan(1);
yout = zeros (length(y0),m+1); % m+1 colonne per farci stare anche la soluzione iniziale
yout(:,1) = y0;
for n = 1:m
    if n == m
        k = tspan(2)-tout(n);
        tout(n+1) = tspan(2);
    else
        tout(n+1) = n*k;
    end
    % Ciclo per costruire le fi
    for i = 1:length(b)
        somma = zeros(length(y0),1);
        for j = 1:i-1
            somma = somma + A(i,j)*f(:,j);
        end
        % lo zero di F(x) sara' f(:,i)
        F = @(x) x - odefun(tout(n)+c(i)*k, yout(:,n)+k*somma + k*A(i,i)*x); 
        J = @(x) eye(length(y0)) - ...
            odeJf(tout(n)+c(i)*k, yout(:,n) + k*somma + k*A(i,i)*x)...
            *k*A(i,i);
        y = yout(:,n); % valore iniziale
        delta = -J(y)\F(y);
        tol =  1e-10;
        while (norm(delta,inf)>tol)
            y = y + delta;
            delta = -J(y)\F(y);
        end
        f(:,i) = y;
    end
    somma = zeros(length(y0),1);
    for j=1:length(b)
        somma = somma + b(j)*f(:,j);
    end
    yout(:,n+1) = yout(:,n) + k*somma;
end
yout = yout';
yrif = yout(end,:);

counter = 0;
for h = 4:1:8
    k = 2^(-h);
    counter = counter + 1;
    m = floor((tspan(2)-tspan(1))/k)+1;
    tout = zeros(1,m+1);
    tout(1) = tspan(1);
    yout = zeros (length(y0),m+1); % m+1 colonne per farci stare anche la soluzione iniziale
    yout(:,1) = y0;
    for n = 1:m
        if n == m
            k = tspan(2)-tout(n);
            tout(n+1) = tspan(2);
        else
            tout(n+1) = n*k;
        end
        % Ciclo per costruire le fi
        for i = 1:length(b)
            somma = zeros(length(y0),1);
            for j = 1:i-1
                somma = somma + A(i,j)*f(:,j);
            end
            % lo zero di F(x) sara' f(:,i)
            F = @(x) x - odefun(tout(n)+c(i)*k, yout(:,n)+k*somma + k*A(i,i)*x);
            J = @(x) eye(length(y0)) - ...
                odeJf(tout(n)+c(i)*k, yout(:,n) + k*somma + k*A(i,i)*x)...
                *k*A(i,i);
            y = yout(:,n); % valore iniziale
            delta = -J(y)\F(y);
            tol =  1e-10;
            while (norm(delta,inf)>tol)
                y = y + delta;
                delta = -J(y)\F(y);
            end
            f(:,i) = y;
        end
        somma = zeros(length(y0),1);
        for j=1:length(b)
            somma = somma + b(j)*f(:,j);
        end
        yout(:,n+1) = yout(:,n) + k*somma;
    end
    yout = yout';
    L(counter) = length(tout);
    Y(counter) = yout(end,:);
    errore(counter) = norm(norm(yrif,inf)-norm(Y(counter),inf),inf);
end

figure
loglog(L,errore,'-*',L,L.^(-1),'k',L,L.^(-2),'r',L,L.^(-3),'g',L,L.^(-4),'y');
legend('errore','pendenza -1','pendenza -2','pendenza -3','pendenza -4');
figure

Questo codice chiaramente ha dei problemi ma non so come modificarlo per farlo funzionare.
Grazie

[malva1]

prova a fare un array lineare e un reshape di counter...così dovrebbe funzionare!