clc;
clear;
format long;

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

% Vetores
x = zeros(N+1,1);
y = 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

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

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