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);

% Condição inicial
x(1) = 0;
y(1) = 1;

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

% Método de Runge-Kutta de Ordem 3
for n = 1:N

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

    k2 = f(x(n) + h/2, y(n) + (h/2)*k1);

    k3 = f(x(n) + h, y(n) - h*k1 + 2*h*k2);

    y(n+1) = y(n) + (h/6)*(k1 + 4*k2 + k3);

    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');