clc;
clear;
format long;

% Dados do problema
a = 0;
b = 2;
N = 10;
h = (b-a)/N;

% Vetores
x = zeros(N+1,1);
y = zeros(N+1,1);
yex = zeros(N+1,1);
erro = zeros(N+1,1);
e
% Condição inicial
x(1) = 0;
y(1) = 1;

% Função
f = @(x,y) x.*y.^(1/3);

% Método de Euler Modificado
for n = 1:N

    k1 = f(x(n),y(n));

    xmeio = x(n) + h/2;
    ymeio = y(n) + (h/2)*k1;

    k2 = f(xmeio,ymeio);

    y(n+1) = y(n) + h*k2;
    x(n+1) = x(n) + h;

end

% Solução exata
yex = ((x.^2 + 2)./2).^(3/2);

% Erro absoluto
erro = abs(y - yex);

% Exibir tabela
fprintf('---------------------------------------------------------------\n');
fprintf('    x         y_aprox         y_exato         erro_abs\n');
fprintf('---------------------------------------------------------------\n');

for i = 1:N+1
    fprintf('%6.2f   %14.10f   %14.10f   %14.10f\n', ...
            x(i), y(i), yex(i), erro(i));
end

fprintf('---------------------------------------------------------------\n');