clc;
clear;
format long;

% Dados do problema
a = 0;
b = 4;
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) 1./(1+4*x.^2) + 0.4*y.^2;

% Método de Runge-Kutta de 4ª ordem
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/2, y(n)+(h/2)*k2);

    k4 = f(x(n)+h, y(n)+h*k3);

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

    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